The Discrete Dantzig Selector: Estimating Sparse Linear Models
via Mixed Integer Linear Optimization

arXiv:1508.01922v3 [stat.ME] 19 Jan 2017

Rahul Mazumder
1

∗1

and Peter Radchenko

†2

Massachusetts Institute of Technology
2
University of Southern California
June, 2016
Abstract

We propose a novel high-dimensional linear regression estimator: the Discrete Dantzig
Selector, which minimizes the number of nonzero regression coefficients subject to a budget
on the maximal absolute correlation between the features and residuals. Motivated by the
significant advances in integer optimization over the past 10-15 years, we present a Mixed
Integer Linear Optimization (MILO) approach to obtain certifiably optimal global solutions to
this nonconvex optimization problem. The current state of algorithmics in integer optimization makes our proposal substantially more computationally attractive than the least squares
subset selection framework based on integer quadratic optimization, recently proposed in [8]
and the continuous nonconvex quadratic optimization framework of [33]. We propose new
discrete first-order methods, which when paired with state-of-the-art MILO solvers, lead to
good solutions for the Discrete Dantzig Selector problem for a given computational budget.
We illustrate that our integrated approach provides globally optimal solutions in significantly
shorter computation times, when compared to off-the-shelf MILO solvers. We demonstrate
both theoretically and empirically that in a wide range of regimes the statistical properties
of the Discrete Dantzig Selector are superior to those of popular `1 -based approaches. We
illustrate that our approach can handle problem instances with p = 10,000 features with certifiable optimality making it a highly scalable combinatorial variable selection approach in
sparse linear modeling.

1

Introduction

We consider the familiar linear regression framework, with response vector y ∈ Rn×1 , model
matrix X ∈ Rn×p , regression coefficients β ∈ Rp×1 and errors  ∈ Rn×1 : y = Xβ + . We assume,
∗

Rahul Mazumder’s research was partially supported by the Office of Naval Research: ONR-024511-00001 and
a grant from the Moore Sloan Foundation. email: rahulmaz@mit.edu
†
Peter Radchenko’s research was partially supported by NSF Grant DMS-1209057. email: radchenk@usc.edu

1

unless otherwise mentioned, that the columns of X, denoted by xj for j = 1, . . . , p, have been
standardized to have zero means and unit `2 -norm. In many modern statistical applications, the
number of variables, p, is larger than the number of observations, n. In such cases, to carry out
statistically meaningful estimation, it is often assumed that the number of nonzero elements in
b which is sparse and serves as a
β is quite small [27]. The task is to obtain a good estimate, β,
good approximation to the underlying true regression coefficient. Of course, the basic problem of
obtaining a sparse model with good data-fidelity is also of interest when the number of observations is comparable to or larger than p. In the sparse high-dimensional setting described above,
two estimation approaches that have been very popular among statisticians and researchers in
related fields are the Lasso [40] and the Dantzig Selector [16]. Both estimators can be expressed as
solutions to convex optimization problems, which can be solved using computationally attractive
procedures [11, 3, 30, 23], and come with strong theoretical guarantees [16, 9, 13]. For reasons
that are explained later in this section, the primary motivation for our investigation in this paper is the Dantzig Selector, which is defined as the solution to the following linear optimization
problem:
min kβk1
s.t.
kX> (y − Xβ)k∞ ≤ δ.
(1)
β

To distinguish this estimator from our proposed approach, we refer to it as the `1 -Dantzig Selector.
This estimator seeks to minimize the `1 -complexity of the coefficient vector, subject to a constraint
on the maximal absolute correlation between the corresponding residual vector and the predictors.
The tuning parameter δ controls the amount of data-fidelity: a small value of δ corresponds to a
good fit, and a larger value of δ leads to heavy shrinkage of the estimated regression coefficients.
[16] point out several reasons as to why the feasibility set in (1) might serve as a good measure
for data-fidelity. In particular, this set is invariant with respect to orthogonal transformations
on the data (y, X). It can also be shown that δ controls1 the residual sum of squares: the
latter can be made arbitrarily close to the minimal least-squares value by decreasing δ. The
`1 -Dantzig Selector, like the Lasso, is used extensively as a model fitting routine to obtain a path
of sparse linear models, as the data-fidelity parameter is allowed to vary [30], and allows a natural
extension to more general response distributions [29]. Note that Problem (1) can be rewritten as
a linear optimization problem and can be solved quite easily for problems with p in the order of
thousands. Under some mild conditions, and even for p much larger than n, the corresponding
estimator achieves a loss within a logarithmic factor of the ideal mean squared error achieved if
the locations of the nonzero coordinates were known [16, 9].
The `1 -Dantzig Selector, however, has limitations. In the presence of highly correlated covariates,
1

More formally, we have:
kX> (y − Xβ)k∞ ≥

λpmin (X) 
1

2np 2

b k22
ky − Xβk22 − ky − Xβ
LS

1
2

,

b is any least-squares solution — see Propowhere, λpmin (X) is the minimum nonzero singular value of X, and β
LS
sition A.1 in [21] for a proof of this result.

2

the estimator tends to choose a dense model, typically bringing in an important variable together
with its correlated cousins, which does not significantly hurt the `1 -norm of the corresponding
coefficient vector. If one increases the data-fidelity threshold δ, the selected model becomes
sparser, however, in the process, important variables might get left out. This is largely due to the
nature of the bias imparted by the `1 -norm, which penalizes both large and small coefficients in
a similar fashion. Similar issues also arise in the case of Lasso [35, 26, 43, 13]. If the `0 -pseudonorm is used instead of the `1 -norm, the aforementioned problems can be ameliorated: given
multiple representations of the model with similar data-fidelity, the `0 -pseudo-norm will always
prefer the most parsimonious representation. In addition, the `0 -pseudo-norm does not shrink
the regression coefficients: once an important variable enters the model, it comes in unshrunk
with its full effect, which, in turn, drains the effect of its correlated cousins and naturally leads
to a sparser model.
Our Proposal. The preceding discussion suggests a natural question: what if we replace kβk1
P
in Problem (1) with kβk0 := pi=1 1(βi 6= 0) – the number of nonzero entries in β? This leads to
the following discrete optimization problem, which also happens to define the estimator that we
propose:
min kβk0
s.t. kX> (y − Xβ)k∞ ≤ δ.
(2)
β

We refer to the above estimator as the Discrete Dantzig Selector. A couple of questions that may
be asked at this point are:
• Is the estimator defined via Problem (2) computationally tractable?
• Does the Discrete Dantzig Selector lead to solutions with superior statistical properties,
when compared to its `1 counterpart?
Addressing these questions and answering them affirmatively is the main focus of this paper.
The objective function in Problem (1), represented by kβk1 , may be thought of as a convexification
of the discrete quantity kβk0 , which counts the number of nonzeros in the regression coefficient
vector β. The corresponding estimator seeks solutions with small `1 -complexity. While this
often leads to sparse solutions, i.e. those with few nonzero coefficients, the sparsity is an indirect
consequence of minimizing kβk1 . The Discrete Dantzig Selector on the other hand, targets
sparsity directly, in its very formulation. Problem (2) can be reformulated as a Mixed Integer
Linear Optimization (MILO) problem — due to the major advances in algorithmic research in
MILO over the past 10-15 years, these methods are widely considered as a mature technology in
a subfield of mathematical programming [41, 28]. Algorithmic advances coupled with hardware
and software improvements have made MILO problems solvable to certifiable optimality for various
problem sizes of practical interest. In this sense, it is perhaps appropriate to perceive MILO as
a computationally tractable tool. The view of computational tractability we adopt here is not
polynomial time tractability, but the ability of a method to provide high quality solutions with

3

provable optimality certificates for problem types that are encountered in practice, in times that
are appropriate for the applications being addressed. Our approach is aligned with an intriguing
recent line of work in computational statistics: the use of Mixed Integer Optimization and, more
broadly, modern optimization techniques to solve certain classes of discrete problems arising in
statistical estimation tasks — see, for example, the recent works of [8, 6]. Further background
on MILO appears in Section 2.1.
In this paper, we bring together recent advances from diverse areas of modern mathematical
optimization methods: first-order techniques in convex optimization and MILO techniques. We
provide a novel unified algorithmic approach that
(a) performs favorably over standalone of-the-shelf MILO solvers applicable for Problem (2), in
terms of obtaining good quality solutions with provable certificates of optimality, and
(b) scales gracefully to problem sizes up to p = 10,000 or even larger.
In an extensive series of experiments with synthetic and real data we demonstrate that our unified
approach solves, to global optimality, instances of Problem (2) with n ≈ 500, p ≈ 100 in seconds,
and underdetermined problems with n ≈ 900, p ≈ 3,000 in minutes. While it takes marginally
longer to provide certificates or guarantees of global optimality, the corresponding times are quite
reasonable: in all the aforementioned instances the certificates of optimality are available within
an hour. Our approach scales to several instances of problems with n ≤ p and p in the range 5,000
to 10,000, delivering optimal solutions in approximately an hour and proving optimality within
at most two days, in all instances. We also find that the statistical properties of our estimates
are substantially better than those of computationally friendlier alternatives, like the `1 -Dantzig
Selector, in terms of both the estimation error and the variable selection properties. Detailed
results appear in Section 7.
Examples. To provide the reader with some intuition, we present a set of three examples, which
illustrate the differences between the solutions to Problems (1) and (2). The following simple
example2 demonstrates how the `1 -based method Dantzig Selector might experience difficulty in
producing a sparse solution in cases where the signal predictors are highly correlated.
Example 1. Let p = n + 1. Take the first feature as x1 = (1, τ, ....., τ )> , take the (i, j)th entry of
the feature matrix, X, as xij = 1 if i = j − 1 and zero otherwise, for i = 1, ..., n and j ≥ 2, and
set y = x1 − x2 .
The `1 -norm of the sparse representation of the response, y = x1 − x2 , equals 2. Note that,
given the available predictors, the response admits only one other exact representation, y =
τ x3 + ... + τ xp . The `1 cost for this dense representation is τ (n − 1), which is lower than the
corresponding value for the sparse representation when τ is small. Consequently, as long as
τ (n − 1) < 2, both the Lasso and the `1 -Dantzig selector select the dense representation of the
2

this example was suggested to us by Emmanuel Candes

4

response. Alternatively, `0 -based methods recover the sparse representation. More specifically,
consider the solution to Problem (2): if the tuning parameter δ is set below τ /(1 + τ ), then the
estimator exactly recovers the sparse representation of the response.
Discrete Dantzig Selector Coefficient Profiles
Example 10
1.0

1

1

1

1

Diabetes data (n=442,p=10)
0.8

1 11

********************************************

1

1

3

3

89
6514317

89
6513217

−4.0

−3.5

−3.0

−2.5

−2.0

2

2

2 2

2

2

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 2* * * * * * * * * * * * * *

2

********************************************

−1.5

−4

−3

−2

0.2

876549131

9

64

4

4

8
1
71

81
61

1
851

2

75
2

6
2
7

2

5

−2.0

−1.5

−1

−3.0

−2.5

9

4
3
9

89
4

0.0

911567 764191358
769131458 8651913147 834

3

39

−0.2

2 22

−1.0

−1.0

2

8765141931

93

−0.4

34567891

3

−0.5

34567891 345678913456789134567891

1
234567891 1234567891

0.0

0.0

4

3

3
9

9

0.4

0.5

0.5

0.6

* * * * * *1* * * * * * * * * * * 1* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

−0.5

Regression Coefficients

1

1.0

Example 1

44
1856 78125 156782

175

1567842

45678912

6
72

6

−1.0

−0.5

********************************************

437619815
12

598417136
12

149386157
21

715648193
12

7143961528

5318472619

96418152731

−0.5
−1.0

8

−0.5

8

********************************************

−4.0

−3.5

−3.0

−2.5

log10 (δ)

−2.0

−1.5

8

8

8

31
79
12541

31
7
421591

13 31
7 7
514912 419125

6

6

6

6

8

8

8

8

6
8

3
1
7
59421

1491527
6

6
8

1
1
7564219
8

0.4

3
14 34
7961512 511872619 9641852731
8

8

3

6
4

0.2

653148917
12

13
957
421

3

9
9
3 3

*************************************************

−4

−3

log10 (δ)

−2

−1

39 3 3 3 3 3 3
9 9 9 9 9 9

46 64 4 4 4
6
8 8 8
8
17 1 1 1 681 18
1 17 17 1 1 15
7
6
7
2 2 2 2 2 72
5

0.0

874631591
21

31 13 13
975 795 759
112 112 12
4 4 4
6
6 6

3

−0.2

715984316
21

13
97
1152
4
6

9 9
3
1

5

−0.4

96415873
21

31
97
1521
4
6

0.0

0.0

75493186
12

0.5

0.5

0.6

*************************************************

−1.0

Regression Coefficients

0.8

1.0

1.0

`1 -Dantzig Selector Coefficient Profiles

4 4

4

1 1 1
185 518 581
6 6 6
2 2 2
7
7 7

4 4
1
581
6
2
7

3 3 3 3
9 9 9 9

3
9 3
9

3
9
4
4
4
3
4
1 1
815 6158 1685 51826 186521 516281 7161528 51847261 964185271
7
6
2 7 7
2
2
7 7
7
4 4

5
5

5

−3.0

−2.5

−2.0

−1.5

−1.0

−0.5

log10 (δ)

Figure 1: Coefficient profiles for the Discrete Dantzig Selector and the `1 -Dantzig Selector, as a function
of the data-fidelity parameter, δ. The dashed vertical lines in the top row of plots indicate the locations
where the number of active variables changes. [Left Panel] corresponds to Example 1; [Middle Panel]
corresponds to Example 10 ; the “true” nonzero coefficients of +1 and −1 are shown as horizontal starred
lines. [Right panel] corresponds to the path for the Diabetes dataset. The numbers overlain on the profiles
indicate the different features.

Figure 1 demonstrates the difference between the Discrete Dantzig Selector and the `1 -Dantzig
selector, by displaying the coefficient profiles for both methods. Note that the profiles for the
Discrete Dantzig Selector are constructed in a piece-wise constant fashion, where for each given
model size, the displayed coefficients are taken from the solution corresponding to the lowest
b ∞ . The left panel of Figure 1 corresponds to Example 1, with
attainable value of kX> (y − Xβ)k
n = 10, p = 11 and τ = 1/2p, which we slightly modified by adding noise to the response: y =
x1 −x2 +. i ’s are independently generated from a centered Gaussian distribution, corresponding
5

to the Signal to Noise Ratio3 (SNR) of 1.3 × 105 .
Now consider Example 10 , which is similar in spirit to Example 1. Here, the first two features,
x1 and x2 , are drawn from a centered bivariate Gaussian distribution with correlation 0.7. The
remaining p − 2 features are drawn from an independent standard Gaussian ensemble. All the
features are standardized to have unit `2 -norm, and the response is generated with SNR=1.4×103 .
The middle panel in Figure 1 displays the corresponding coefficient profiles with n = 10, p = 12.
The Discrete Dantzig Selector exactly recovers the true model for a wide range of the tuning
parameter, δ. As δ is decreased, and noise variables come into the model, their coefficients
remain highly shrunk, while the coefficients for the signal variables remain near their true values.
On the other hand, the `1 -Dantzig selector is unable to recover the true model, and produces a
large estimation error for all values of the tuning parameter.
The right panel in Figure 1 corresponds to the well known Diabetes dataset [19], where n = 442
and p = 10; here all the variables including the response were standardized to have unit `2 -norm
and zero mean. Note that, despite some similarities, the sequences of predictors entering the
model are different for the two approaches.
We note that instead of formulation (2), one may prefer to minimize the data-fidelity term subject
to a constraint on the number of nonzeros in β:
min kX> (y − Xβ)k∞
β

s.t.

kβk0 ≤ k.

(3)

The framework developed in this paper may be adapted to Problem (3). In this paper, however,
we focus on Problem (2).
Context and Related Work. A primary motivation of our work is derived from the recent
work on the least squares variable selection problem [8], where the authors study
min
β

1
ky − Xβk22
2

s.t.

kβk0 ≤ k,

(4)

using mixed integer convex quadratic optimization (MIQO) methods. While certain statistical
properties of solutions from Problems (2) and (4) are comparable; we observed in our computational experiments (see Section 7.1) that Problem (2) is orders of magnitude faster (often by
a factor of hundreds) than Problem (4) in obtaining solutions with certificates of global optimality. In addition, a MILO formulation for Problem (2) consumes much less memory than a
comparable MIQO formulation for Problem (4). The aforementioned computational superiority
of MILO over MIQO should not come as a surprise. Indeed, it is quite well known in the integer
programming community (see, for example, the nice review papers [28, 14]) that current algorithms for MILO problems are a much more mature technology than MIQO. Recently, [33] proposed
MIPGO for nonconvex penalized least squares regression based on purely continuous nonconvex
3

For a model generated as yi = µi + i , i = 1, . . . , n; we define SNR as follows: SNR= Var(µ)/Var().

6

quadratic optimization. We demonstrate the substantial superiority of our proposal over MIPGO
(in Section 7.1) in obtaining high quality statistical solutions within a given computational budget.
Thus, the superior computational scalability of the corresponding optimization methods forms a
principal motivation to study Problem (2), as an effective estimation procedure for sparse linear
regression. In addition, from a statistical viewpoint, in terms of estimating sparse regression
models subject to good data-fidelity, the Discrete Dantzig Selector may be perceived as a natural,
interpretable and useful alternative to least squares with variable selection, Problem (4) — in the
same way as the `1 -Dantzig Selector may be viewed as an appealing alternative to Lasso.
Contributions. Our main contributions can be summarized as follows:
1. We propose a new high-dimensional linear regression estimator: the Discrete Dantzig Selector, which minimizes the number of nonzero regression coefficients, subject to a budget on
the maximal absolute correlation between the features and the residuals. We show that the
estimator can be expressed as the solution to a MILO problem, a computationally tractable
framework that delivers certifiably optimal global solutions; and is computationally more
scalable than the recently proposed methods in [8, 5, 33].
2. We develop new discrete first-order methods, motivated by recent algorithmic developments
in first-order continuous convex optimization, to obtain high quality feasible solutions for
the Discrete Dantzig Selector problem. These solutions are passed onto MILO solvers as
warm-starts. Our proposal leads to advantages over the off-the-shelf state-of-the-art integer programming algorithms in terms of (a) obtaining superior upper bounds for a given
computational budget and (b) aiding MILO solvers in obtaining tighter lower bounds and
hence improved certificates of optimality. Exploiting problem specific information, we also
propose enhanced MILO formulations, which further improve the algorithmic performance
of MILO solvers.
3. We characterize the statistical properties of the Discrete Dantzig Selector and demonstrate
both theoretically and empirically its advantages over `1 -based approaches. Our results also
apply to approximate solutions for the Discrete Dantzig Selector optimization problem.
4. Our approach obtains optimal solutions for p ≈ 500 in a few minutes, p ≈ 3,000 within
fifteen minutes and for problems with p = 10,000 in an hour. Certificates of optimality
are obtained at the expense of higher computation times—for instances with p ≈ 500 they
are obtained within half-hour, for p ≈ 3,000 they are achieved around an hour and for
p = 10,000 the certificates arrive in the range from three to forty hours. To the best of our
knowledge, we present herein, the largest problem instances in subset selection for which
certifiably optimal solutions can be obtained.
Roadmap. The remainder of the paper is organized as follows. Section 2 describes the op-

7

timization methodology behind the proposed approach, and discusses its connections with the
`1 -Dantzig Selector optimization problem. In Section 3, the statistical properties of the Discrete Dantzig Selector are analyzed from a theoretical point of view; the results are compared
to the `1 -Dantzig Selector and the Lasso. The framework of discrete first-order methods is described in Section 4. Additional discussion of MILO formulations together with problem specific
enhancements, is presented in Section 5. Section 6 gathers numerical results on the computational
performance of our algorithms in a variety of settings. An empirical analysis of the statistical
properties of the Discrete Dantzig Selector is conducted in Section 7. Some technical details are
provided in the Appendix.

2

Overview of the Proposed Methodology

Herein, we introduce and summarize the general aspects of the proposed methodology. Further
details and enhancements are provided in Sections 4 and 5.

2.1

Mixed Integer Linear Optimization (MILO) Preliminaries

The general form of a MILO problem is as follows:
min

a> α

s.t.

Aα ≤ b

α

αi ∈ {0, 1},

i = 1, . . . , m1

αj ≥ 0,

j = m1 + 1, . . . , m,

where a ∈ Rm×1 , A ∈ Rd×m and b ∈ Rd×1 are the problem data, the symbol “≤” denotes
element-wise inequalities, and we optimize over (α1 , . . . , αm ) := α ∈ Rm containing both discrete
(αi , i = 1, . . . , m1 ) and continuous (αi , i = m1 + 1, . . . , m) variables. For background on MILO,
we refer the reader to [7, 31]. Some modern integer optimization solvers include Cplex, Glpk,
Gurobi, Knitro, Mosek, Scip — see also [32].
As already alluded to in Section 1, there has been significant progress in the theory and practice
of MILO over the past fifteen to twenty years. Specifically, the computational power of MILO
solvers has undergone impressive advances over the past twenty-five years — the cumulative
machine-independent speedup factor in MILO solvers between 1991 and 2015 is estimated to be
780,000 [10]. This progress can be attributed to the inclusion of both theoretical and practical
advances into MILO solvers. Some of the main factors responsible for this speedup are advances
in cutting plane theory, improved heuristic methods, disjunctive programming for branching
rules, techniques for preprocessing MILOs, using linear optimization as a black box to be called
8

50
40

●

●

100

Upper Bounds
Lower Bounds
Global Minimum

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

20

30

Bounds

●

●
●
●

10

●

20
0

0

●
●

500

1000

1500

●
●

●
●
●
●
●

0

●
●
●
●

●

●
●
●

●
●
●
●

80

60

●

60

●

40

MILO Optimality Gap (in %)

b 0 = 41)
Diabetes Dataset (n = 442, p = 64, kβk

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

0

Time (secs)

500

1000

1500

Time (secs)

Figure 2: Typical evolution of a MILO algorithm for Problem (2), as a function of time. [Left Panel]
displays the progress of Upper Bounds (UB) and Lower Bounds (LB) for the optimal value of the objective
function. The upper bounds, which correspond to feasible solutions for Problem (2) are seen to stabilize at
the global minimum within a few seconds. The lower bounds provide a certificate of how far the current
solution might be from the global solution — these bounds progressively improve as the MILO algorithm
explores more nodes in the branch and bound tree. Observe that the certificate of global optimality arrives
at a later stage, even though the algorithm finds the global solution very quickly. [Right Panel] displays
the evolution of the corresponding MILO Optimality Gap (in %), defined as (UB - LB)/UB, with time.

by MILO solvers, and improved linear optimization methods [10]. In addition, there have been
substantial improvements in hardware speed: the overall hardware speedup from 1993 to 2015 is
approximately estimated to be 105.75 ∼ 570,000 [1]. When both hardware and software advances
are combined, the overall speedup for MILO problems is estimated to be around 450 billion! One
attractive feature of MILO solvers, which is a stark contrast to heuristic approaches, is that the
former provide (a) feasible solutions, which are also upper bounds to the minimum objective value
and (b) lower bounds for the optimal value of the objective function. As a MILO solver makes
its way to the global optimum, the lower bounds become tighter, thereby providing improved
certificates of sub-optimality (see Figure 2 for an illustration). This aspect of MILO solvers is
quite useful, especially if one decides to stop the solver before reaching the global optimum. In
the modern day world, MILO plays a key role in various impactful application areas of operations
research: revenue management, air-traffic control, scheduling and matching tasks, production
planning and others [42, 7]. In this paper, we show how the power of MILO can be used in
the context of a problem of fundamental importance in statistics, namely, sparse linear model
estimation — we build upon recent line of work in computational statistics, at the interface of
modern discrete optimization and fundamental techniques in statistical modeling [6, 8].

9

2.2

MILO formulations for the Discrete Dantzig Selector

Assuming without loss of generality that Problem (2) has a minimizer which is bounded, it can
be obtained by solving
kβk0

Γ1 := min
β

kX> (y − Xβ)k∞ ≤ δ

s.t.

(5)

kβk∞ ≤ MU ,
where, MU is a large but finite number [7]. We present a MILO formulation for Problem (5) (and
also, Problem (2))
p
P
zi
min
β,z

i=1

s.t.

−δ ≤ dj − hqj , βi ≤ δ,

j = 1, . . . , p

−MU zj ≤ βj ≤ MU zj ,

j = 1, . . . , p

zj ∈ {0, 1},

j = 1, . . . , p,

(6)

where the optimization variables are z (binary) and β (continuous); the problem data consists of
d := X> y = (d1 , . . . , dp )> and Qp×p := X> X = [q1 , . . . , qp ]. Formulation (6) is often referred
to as a “Big-M” formulation due to the presence of the parameter MU . The binary variable zi
controls whether βi is zero or not: if zi = 0 then βi = 0 and if zi = 1 then βi is free to vary in
P
the interval [−MU , MU ]. The objective function pi=1 zi controls the number of nonzeros in the
model. Figure 2 shows the performance of the above MILO formulation on the Diabetes dataset [19]
with n = 442, p = 64 (here, we mean-centered and scaled y, xi ’s to have unit `2 -norm).
Formulation (6) has intriguing connections to the `1 -Dantzig Selector: the binary variables zi ∈
{0, 1} in Problem (6) can be relaxed into continuous variables zi ∈ [0, 1], leading to:
Γ2 := min
β

s.t.

1
MU

p
P

|βj |

i=1

(7)

kX> (y − Xβ)k∞ ≤ δ
−MU ≤ βj ≤ MU ,

j = 1, . . . , p.

Problem (7) modifies the `1 -Dantzig Selector problem:
Γ3 := min

1
kβk1
MU

s.t.

kX> (y − Xβ)k∞ ≤ δ,

by including an `∞ -constraint on β. It follows from the above that: Γ3 ≤ Γ2 ≤ Γ1 . The last
inequality is typically strict, and, depending upon the data, the gap between the values, as well as
between the corresponding optimal solutions, can be substantial, as illustrated in Figure 1. The
above discussion provides another viewpoint for explaining the differences between the `1 -Dantzig

10

Selector and the Discrete Dantzig Selector estimators. If MU is taken to be large enough then
Γ3 = Γ2 .
We note that the constraint kβk0 ≤ k can also be expressed using Specially Ordered Sets [7] as
an alternative to the “Big-M” formulation (6). This can be used when the user does not wish to
specify a-priori any bound MU on the regression coefficients and/or the coefficients have widely
varying amplitudes. We discuss this further in Section 5.
We emphasize that MU appearing in formulation (6) should not be interpreted as a statistical
tuning parameter — it appears from a purely algorithmic viewpoint and, as we saw, has interesting
connections to the `1 -Dantzig Selector problem. MU might be taken to be arbitrarily large (but
finite) to obtain a solution to Problem (2). In Section 5 we describe several data driven methods
to estimate MU , and we also discuss other structured formulations of Problem (6), which lead to
improved algorithmic performance: they deliver tighter computational lower bounds in smaller
amounts of time. We now proceed towards an analysis of the statistical properties of the Discrete
Dantzig Selector and investigate its comparative advantages over its `1 -counterpart.

3

Statistical Properties: Theory

In this section we study the statistical properties of the Discrete Dantzig estimator. In particular, we characterize its connections with the Best Subset selection estimator in the case of the
orthonormal design, we investigate its oracle properties in the classical asymptotic regime and,
finally, we analyze its estimation, prediction and variable selection performance in the high dimensional setting. To improve readability, all the technical proofs are presented in Section A in
the Appendix.

3.1

Orthonormal Design

Here we assume x>
j xk equals 1 if j = k and 0 otherwise. Note that such an assumption requires
b which solves
n ≥ p. Our goal is to compare and connect the Discrete Dantzig estimator, β,
BS
b , which solves (4).
the optimization problem (2), with the Best Subset selection estimator, β
In particular, we want to understand the relationship between the tuning parameter δ in the
Discrete Dantzig optimization problem and the tuning parameter k, which controls the `0 norm
of the Best Subset solution. Note that Problem (2) does not generally have a unique optimizer,
b to refer to just one of the Discrete Dantzig solutions.
so we use β
Define cj = x>
j y for j = 1, ..., p. To simplify the presentation, and without loss of generality,
suppose that the predictors are indexed in such a way that |c1 | ≥ |c2 | ≥ ... ≥ |cp |. Suppose also
that ck 6= ck+1 , to ensure uniqueness of the Best Subset solution.

11

Theorem 1.

1. The Best Subset selection estimator is uniquely defined as
b BS = (c1 , . . . , ck , 0, . . . , 0).
β

2. Suppose that δ ∈ [ck+1 , ck ). Then, the set of all the Discrete Dantzig solutions is
{(c1 + u1 , ..., ck + uk , 0, ..., 0) , |uj | ≤ δ, j = 1, ..., k} .
It follows that each Best Subset estimator is a Discrete Dantzig solution for an appropriately
chosen δ. The coefficients of both estimators are obtained by the hard thresholding of the covariances cj , however, the nonzero coefficients of the Discrete Dantzig estimator are allowed to
deviate from the value of the covariance by an amount bounded above by δ.

3.2

Fixed p asymptotics

To avoid confusion, we refer to the true coefficient vector as β ∗ . In this subsection we treat the
number of predictors, p, as fixed and let the number of observations, n, tend to infinity. The
standard assumption for deriving asymptotic results in this setting is that β ∗ does not depend
on n, and n−1 X> X converges to a non-singular covariance matrix C. We rewrite the above
assumption to be consistent with the scaling kxj k2 = 1 for j = 1, ..., p, which is used throughout
∗
this paper. Thus, we require that X> X converges to C as n tends to infinity, and β ∗ = n1/2 β̃ ,
∗
for some fixed vector β̃ . We also impose a standard assumption that i are i.i.d. with zero mean
and finite variance. Given an index set J ⊆ {1, . . . , p} we write XJ for the sub-matrix of X that
consists of the columns identified by J. Let J ∗ denote the support of the vector β ∗ . We define
b O , as the least-squares estimator computed using only the true predictors.
the Oracle estimator, β
b O equals J ∗ , and β
b O∗ = (X>∗ XJ ∗ )−1 X>∗ y.
In other words, the support of β
J
J
J
Theorem 2. Let δ → ∞ and δ = o(n1/2 ). Suppose that matrix C is invertible. Then, with
probability tending to one,
1. The support of each solution to the Discrete Dantzig optimization problem equals J ∗ ;
b O , belong to the set of
2. Both the true coefficient vector, β ∗ , and the Oracle estimator, β
solutions to the Discrete Dantzig optimization problem.
Consider the polished version of the Discrete Danzig estimator, which is defined as follows: given
b P , is set
b the support of the polished estimator, β
a Discrete Dantzig solution with support J,
b P is defined as the least-squares estimator using the corresponding
b on its support β
equal to J;
b Pb = (X> X b)−1 X> y. Note that the value of β
b P generally depends on the
predictors, i.e. β
J
Jb J
Jb
the choice of the Discrete Dantzig solution. However, this choice becomes irrelevant under the
setting of Theorem 2. More specifically, with probability tending to one, every polished estimator
coincides with the Oracle estimator.
12

bP = β
b O holds with probability
Corollary 1. Under the assumptions of Theorem 2, equality β
tending to one.
b P satisfies the oracle property in the sense of Fan and Li [20].
Consequently, estimator β

3.3

High Dimensional Setting

Here we focus on the case where p is large, possibly much larger than n. We discuss the properties of the global as well as approximate solutions to Problem (2). We also comment on the
estimator obtained from the closely related Problem (3), in which kβk0 is constrained, rather
than minimized. We assume that the error terms in the underlying linear model are mean zero
Gaussian with variance σ 2 and, as before, use β ∗ to refer to the true regression coefficient vector.
We start with some notation. For every vector θ ∈ Rp and index set J ⊆ {1, ..., p} we write θ J
for the sub-vector of θ determined by J.
Definition 1. Given positive integers k and m, such that m ∈ [k, p − k], and a positive c0 , let
γ(k) =

min

θ6=0, kθk0 ≤k

kXθk2
kθk2

(
κ(k, c0 ) = min

|J0 |≤k

min

θ6=0, kθ J c k1 ≤c0 kθ J0 k1
0

(
κ(k, c0 , m) = min

|J0 |≤k

min

θ6=0, kθ J c k1 ≤c0 kθ J0 k1
0

kXθk2
kθ J0 k2

)

kXθk2
kθ J01 k2

)
,

where J0 ⊆ {1, ..., p} and J01 := J0 ∪ J1 , with J1 identifying the m largest (in magnitude) coordinates of θ outside of J0 .
We use s∗ to denote kβ ∗ k0 . As we discuss in the next subsection, quantities [κ(s∗ , c0 )]−1 and
[κ(s∗ , c0 , m)]−1 , for m ≥ s∗ and c0 ≥ 1, appear in the error bounds for the Lasso and the original
Dantzig selector, while [γ(2s∗ )]−1 appears in the bounds for Discrete Dantzig Selector. The
following result establishes some useful relationships for these quantities.
Proposition 1. For all positive integers k and m, with m ∈ [k, p−k], and all c0 ≥ 1 the following
√
holds: γ(2k) ≥ κ(k, c0 )/ 2 and γ(2k) ≥ κ(k, c0 , m).
Recall the setting of Example 1. When τ (n − 1) < 2, the `1 methods, such as the original Dantzig
selector, fail to recover the sparse representation of the response. Note also that κ(2, c0 ) =
κ(2, c0 , m) = 0, for m ≥ 2 and c0 ≥ 1. On the other hand, γ(4) > 0 for p > 4, and the Discrete
Dantzig Selector succeeds in recovering the correct sparse representation, for every sufficiently
small value of the tuning parameter, δ.
The following theorem establishes several useful bounds for the Discrete Dantzig Selector.
p
b solves optimization problem (2) for δ = σ 2(1 + a) log p, where a ≥
Theorem 3. Suppose that β

13

√
0. The following bounds hold with probability bounded below by 1 − (pa π log p)−1 :
b 0 ≤ s∗
kβk
b − β ∗ k1 ≤ 4 [γ(2s∗ )]−2 s∗ δ
kβ
b − β ∗ k2 ≤ 8 [γ(2s∗ )]−4 s∗ δ 2
kβ
2
b − β ∗ )k2 ≤ 8[γ(2s∗ )]−2 s∗ δ 2 .
n−1 kX(β
2
Remark.

It follows from the proof of Theorem 3 that the above result

(i) holds uniformly over the set {β ∗ : kβ ∗ k0 ≤ s∗ };
(ii) also holds for the solution to Problem (3) with k = s∗ .

We now compare the above bounds to those established for the popular `1 -based approaches.
Under the assumed scaling of the predictors, and for every positive integer m, such that m ∈
[s∗ , p − s∗ ], Theorem 7.1 in [9] gives the following error bounds for the `1 -Dantzig Selector estib :
mator, β
DS
b − β ∗ k1 ≤ 8 [κ(s∗ , 1)]−2 s∗ δ
kβ
DS

2
p
b
kβ DS − β ∗ k22 ≤ 16 1 + s/m [κ(s∗ , 1, m)]−4 s∗ δ 2

(8)

b − β ∗ )k2 ≤ 16[κ(s∗ , 1)]−2 s∗ δ 2 .
n−1 kX(β
DS
2
By Proposition 1, the right hand sides of the above inequalities are at least as large as the
corresponding bounds in Theorem 3. Moreover, the differences in the two sets of bounds can
potentially be quite significant. Consider the setting of Example 1 for illustration. The upper
bounds in Theorem 3 are finite for p > 4, while the three bounds in display (8) are infinite.
Examining Theorem 7.2 in [9], we conclude that the corresponding error bounds for the Lasso
are at least as large as those given in display (8). The same result also provides an upper bound
b
on the `0 -pseudo-norm of the Lasso estimator, β
Lasso :
b
kβ
Lasso k0 ≤

64φmax ∗
s ,
[κ(s∗ , 3)]2

where φmax is the maximum eigenvalue of the matrix X> X. Note that the right-hand side of
the above bound is infinite in the setting of Example 1. In general, this upper bound is at
least 64 times as large as the one for the Discrete Dantzig Selector estimator. We informally
summarize the above findings as follows: when compared to the `1 -based approaches, the Discrete Dantzig Selector satisfies as good or better estimation and prediction error bounds, while
achieving significantly higher level of sparsity.
We can sharpen the bounds in Theorem 3 by making them dependent on the support of β ∗ . More

14

specifically, given an index set J ∗ we define
γ̃(J ∗ ) =

min

J⊂{1,...,p},|J|=2s∗ ,J⊇J ∗

kXθk2
.
θ6=0, θ J c =0 kθk2
min

Then, Theorem 3 holds with γ(2s∗ ) replaced by γ̃({k : βk∗ 6= 0}), and the corresponding result is
uniform over β ∗ .
The following corollary to Theorem 3 shows that the Discrete Dantzig Selector successfully recovers the support of the true coefficient vector, provided the nonzero coefficients are appropriately
bounded away from zero. Define |β ∗ |min = min{|βk∗ |, βk∗ 6= 0}.
p
Corollary 2. If |β ∗ |min > 4σ [γ(2s∗ )]−2 (1 + a)s∗ log p, then the estimator from Theorem 3
√
exactly recovers the support of β ∗ , with probability bounded below by 1 − (pa π log p)−1 .
b that is a feasible solution to the optimization problem (2), but
We now consider an estimator β
b together with a lower
not necessarily the optimal solution. Recall that our algorithms produce β
bound on the minimum value of the objective function, kβk0 . We denote this lower bound by
b 0 is within a prespecified
sbLB . The next result shows that if the algorithm is stopped when kβk
multiplicative factor of sbLB , the bounds from Theorem 3 continue to hold after an appropriate
adjustment. The corresponding proof follows the argument in the proof of Theorem 3, making
only minor modifications.
b is a feasible solution to the optimization problem (2), corresponding
Theorem 4. Suppose that β
p
b 0 ≤ (1 + ψ)b
sLB . Then, the following bounds
to δ = σ 2(1 + a) log p, where a ≥ 0, such that kβk
√
a
−1
hold with probability bounded below by 1 − (p π log p) :
b 0 ≤ (1 + ψ)s∗
kβk
b − β ∗ k1 ≤ 2(2 + ψ) [γ ([2 + ψ]s∗ )]−2 s∗ δ
kβ
b − β ∗ k2 ≤ 4(2 + ψ) [γ([2 + ψ]s∗ )]−4 s∗ δ 2
kβ
2
b − β ∗ )k2 ≤ 4(2 + ψ)[γ([2 + ψ]s∗ )]−2 s∗ δ 2 .
n−1 kX(β
2
Note that the constant ψ is typically quite small in practice, for example, 0.1 (see the right panel
in Figure 2 for an illustration of the evolution of ψ over time for the Diabetes dataset.) Thus,
the corresponding effect on the error bounds is generally minor.

4

Obtaining Good Solutions via Discrete First-Order Methods

In this section we propose new algorithms, referred to as discrete first-order methods, which deliver
good upper bounds for Problem (2). It is important to note that unlike the MILO framework,
these algorithms do not provide lower bounds. Instead, the solutions obtained by our methods are
passed to MILO solvers as warm-starts. The proposed algorithms are inspired by recent advances
15

of first-order methods in convex optimization [37, 36, 38], and can be viewed as their nonconvex
adaptations. We summarize their key advantages:
• They provide excellent upper bounds to Problem (2) with low computational cost, time and
memory requirements.
• MILO solvers accept these solutions as warm-starts and consequently improve upon them.
This hybrid approach outperforms the stand-alone capabilities of an off-the-shelf MILO solver,
producing high quality upper bounds in amounts of time that are orders of magnitude
smaller.
• The solutions obtained can be used to improve the overall run-time of MILO solvers, including
certificates of global optimality.
We validate our proposed methods on several synthetic and real-data datasets.

4.1

Discrete First-Order Methods

Problem (2) involves the minimization of a discontinuous objective function over a polyhedral
set. Thus, it is not directly amenable to simple proximal gradient type algorithms [36, 37, 8].
We propose two algorithms: Algorithm 1 (see Section 4.1.1) and Algorithm 2 (see Section 4.1.2),
both of which can be used as stand-alone solvers for obtaining good quality upper bounds to
Problem (2). We also present a hybrid method, Algorithm 3, which combines the strengths of
both Algorithms 1 and 2, by using the solution obtained from Algorithm 1 as an initialization to
Algorithm 2. In our experiments Algorithm 3 showed the best empirical performance. Section 6
presents numerical results illustrating the performance of our framework. We emphasize that
Algorithms 1—3 only provide good upper bounds, they do not certify the quality of solutions via
lower bounds. A main purpose of these algorithms is to provide good quality upper bounds to
initialize the MILO solvers — the latter, in turn, are often found to improve upon the upper bounds
obtained from these first-order algorithms, at the cost of more (but still reasonable) computation
times.

4.1.1

The Variable Splitting Method

We present our first discrete first-order method based on a classical method in nonlinear optimization: the Alternating Direction Method of Multipliers (aka ADMM) [4] popularly used in
the context of convex optimization–we refer the reader to [12] for a nice exposition on this topic.
We choose this method because of its simplicity and good performance in practice, as seen in our
experiments. To apply this algorithm, we decouple the feasible set {β : kX> (y−Xβ)k∞ ≤ δ} and
the discontinuous function β 7→ kβk0 . Observe that Problem (2) can be equivalently rewritten

16

as:
min

kβk0

s.t.

kX> (y − Xα)k∞ ≤ δ

α,β

(9)

α = β.
We consider the Augmented Lagrangian given by:
Lλ (β, α; ν) := kβk0 +

λ
kβ − αk22 + hν, α − βi
2

(10)

for some value of λ > 0, where, ν may be thought of as a “dual” variable4 that along with λ
controls the proximity between α and β. The ADMM procedure leads to the following update
sequence:
β k+1 ∈ arg min Lλ (β, αk ; ν k )

(11)

αk+1 ∈

(12)

β

arg min
α:kX> (y−Xα)k∞ ≤δ

Lλ (β k+1 , α; ν k )


ν k+1 =ν k + λ αk+1 − β k+1 .

(13)

Step (11) can be performed via a hard thresholding operation [18] and (12) involves a simple
projection onto the polyhedron {α : kX> (y − Xα)k∞ ≤ δ}, which can be done efficiently, as
detailed in Section B.1 (Appendix).
Algorithm 1
(1.) Input (β 1 , α1 , ν 1 ), choose λ > 0, and repeat the following steps until convergence.
(2.) Update (β k , αk , ν k ) to (β k+1 , αk+1 , ν k+1 ) via (11)–(13).
(3.) Stop if kβ k+1 − β k k2 ≤ τ1 kβ k k2 and5 kβ k − αk k2 ≤ τ2 max{kβ k k2 , kαk k2 }, otherwise go to
Step 2.
We found Algorithm 1 to work quite well in our experiments. The algorithm may be sensitive to
the choice of λ — affecting the solution and the time until convergence. We recommend using
multiple values of λ, and choosing the best solution among them. In Section 4.1.3 we address
these shortcomings and describe modifications that lead to improvements in practice.
4
Following the terminology in [12], if instead of kβk0 we had a convex function, then ν would be a dual variable,
and its corresponding update step (13) would be the dual update. We will with a slight abuse of terminology use
the term “dual” here.
5
Here, τ1 , τ2 are tolerances for convergence, typically, taken to be equal and set to 10−4 .

17

4.1.2

Sequential Linear Optimization

We now describe another nonlinear optimization algorithm for obtaining upper bounds for Problem (2), motivated by ideas popularly used in nonconvex penalized regression (see, for example,
[35] and references therein). Let us consider a family of nonconvex functions, ργ (|β|), parametrized
by γ ∈ (0, γ̄], such that γ = γ̄ corresponds to ργ (|β|) = |β|, and, as γ decreases to 0, ργ (|β|)
becomes a progressively better approximation to 1(β 6= 0). In other words,
p
X

kβk0 = lim

γ→0+

ργ (|βi |).

(14)

i=1

We make the following assumption about ργ (·):
Assumption (A): ργ (β) is symmetric in β around zero and continuous. Let ργ (β) ≥ 0 and
ργ (0) = 0. For every γ, the map |β| 7→ ργ (|β|) is concave and differentiable on (0, ∞).
Some popular choices of ργ (·) are ργ (t) = log( γt + 1)/ log( γ1 + 1) and ργ (t) = tγ for t ≥ 0. We refer
the reader to [35] (and references therein) for more context and examples of nonconvex penalty
functions used in sparse linear regression.
We propose to compute good upper bounds for the following continuous nonconvex optimization
problem:
p
X
min
h(β) :=
ργ (|βi |)
β
(15)
i=1
s.t.

kX> (y − Xβ)k∞ ≤ δ,

especially for γ ≈ 0+. In light of (14), this leads to good upper bounds for Problem (2). Note
e
that the concavity of |β| 7→ ργ (|β|) leads to the following upper bound (for all β and β):
h(β) =
≤

p
X
i=1
p
X

|i=1

ργ (|βi |)
ργ (|βei |) +

p
X



ρ0γ (|βei |) |βi | − |βei | ,

(16)

i=1

{z

}

e
:=h̄(β;β)

where ρ0γ (·) denotes the derivative of |β| 7→ ργ (|β|), with the convention that ρ0γ (0) = ∞ if the
derivative is unbounded as |β| → 0+. Inequality (16) suggests that we sequentially minimize an

18

upper bound to the objective function in (15). This leads to the following iterative scheme:
p
P

β k+1 ∈ arg min
β

s.t.

i=1

kX> (y

ρ0γ (|βik |)|βi |

(17)

− Xβ)k∞ ≤ δ,

where we assume, without loss of generality, that β 1 is feasible for Problem (17). The above
sequential approximation of the function h(β) is similar to the popular reweighted `1 -minimization
method, used in signal processing [15] for sparse linear model estimation.
We now present a simple finite time convergence rate of the iterative process (17) in terms of
reaching an approximate first order stationary point of Problem (15). Towards this end, we
introduce the following quantity:
∆(θ) := min
β

s.t.

p
P

ρ0γ (|θi |) (|βi | − |θi |)

i=1

kX> (y

(18)

− Xβ)k∞ ≤ δ,

which we use to define a first-order stationary point for Problem (15).
b is feasible for Problem (15). We say that θ
b is a first-order stationary
Definition 2. Suppose θ
b = 0. For φ > 0, θ
b is said to be an φ accurate first-order stationary
point for Problem (15), if ∆(θ)
b ≥ −φ.
point if ∆(θ)
b is
We note that ∆(θ) (assuming that θ is feasible for Problem (15)) is a measure of how far θ
from a first order stationary point of Problem (15)–if ∆(θ) < 0, then the current estimate θ can
be improved, if ∆(θ) = 0 then the solution cannot be improved via update (17). We refer the
reader to Section B.4 for a more detailed explanation. The following theorem (for a proof see
Section B.4) presents a finite time convergence rate analysis of the sequence (17) to a first-order
stationary point for Problem (15).
Theorem 5. Consider Problem (15) for a fixed γ > 0, with Assumption (A) on ργ (·) in place.
The update sequence β k , defined via (17), leads to a decreasing sequence of objective values for
Problem (15): h(β k+1 ) ≤ h(β k ) for all k ≥ 1. In addition, for every K > 0 we have the following
finite-time convergence rate:

n
o
1 
h(β 1 ) − b
h ,
−∆(β k ) ≤
1≤k≤K
K
min

where the sequence of objective function values satisfies h(β k ) ↓ b
h as k → ∞.
We emphasize that the result in Theorem 5 pertains to the performance of the sequence (17) as
a numerical optimization scheme, and has no direct implication on the statistical properties of
the sequence. Theorem 5 implies that for any φ > 0, it takes at most K = O( φ1 ) many iterations
to reach a φ-accurate first-order stationary point, i.e., there exists a 1 ≤ k ∗ ≤ K such that
19

∗
b , an upper bound for Problem (15) for a
∆(β k ) > −φ. The sequence β k leads to an estimate β
γ
fixed γ. Since our intent is to obtain a good solution to Problem (2), we make use of property (14).
This suggests that we obtain a good upper bound to Problem (15) for a small value of γ ≈ 0+.
Instead of applying iteration (17) for a pre-specified (small) value of γ, we recommend using a
continuation strategy in practice. We take a sequence of decreasing values of γ ∈ {γ1 , . . . , γN },
b as a warm-start for obtaining a good solution (upper bound) to
where γi > γi+1 . We use β
γi
Problem (15) for a smaller value of γ = γi+1 . In our numerical experiments, this continuation
strategy seems to work well. The method is summarized below.

Algorithm 2
(1.) Take a decreasing sequence of γ values {γ1 , . . . , γN }; initialize with β 0 = 0; and fix a value
of Tol = 10−5 (say). Set κ = 1 and γ = γκ .

(2.) Use the update sequence rule (17) until some convergence criterion is met: −∆(β k ) < Tol.
b denote the estimate of β, upon convergence.
Let β
γ
b . If κ ≤ N , then goto Step 2. If κ > N , exit with β
b
(3.) Set κ ← κ + 1, γ = γκ and β 0 ← β
γ
γ
as an upper bound to Problem (2).
The linear optimization Problem (17) can be solved quite efficiently using simplex methods. For
larger problems, i.e. p larger than a few thousand, we recommend using modern first-order
method as described in Section B.2. Since Algorithm 2 requires solving several instances of
related problems of the form (17), the warm-start capabilities of simplex methods and first-order
methods lead to computational benefits.

4.1.3

Algorithm 3: Combining the Strengths of Algorithm 1 and Algorithm 2

In our empirical studies we observed that Algorithm 1 is more effective in obtaining good upper
bounds than Algorithm 2 for a given time limit. Algorithm 2, on the other hand, has stronger
convergence guarantees than Algorithm 1. Algorithm 1 leads to an estimate of β that is sparse
but approximately satisfies6 the feasibility constraint of Problem (2). Algorithm 2, in contrast,
leads to solutions that are both sparse and feasible — these advantages make Algorithm 2 an
important tool in our framework. We propose to combine the best features of Algorithms 1 and
2 to develop a hybrid variant: Algorithm 3, which we recommend to use in practice. Algorithm 3
b (1) , to create a
is simple but very effective: it uses the solution obtained from Algorithm 1, say, β
b (1) , and then applies Algorithm 2 on this set
set I ⊂ {1, . . . , p}, which includes the nonzeros in β
I — the details of this method are presented in Section B.5 (in the Appendix).
6

This is because Algorithm 1 delivers a pair, α, β, which are approximately equal: α ≈ β; α is feasible for
Problem (2) but need not be exactly sparse; β, on the other hand, is sparse but approximately feasible.

20

5

Structured MILO Formulations and Certificates of Optimality

This section is dedicated to enhancements of the basic Discrete Dantzig Selector formulation (6),
presented in Section 2.2. These are particularly useful in delivering tighter lower bounds, thereby
providing certificates of global optimality in shorter times.
Note that formulation (6) requires the specification of MU large enough to include the solution
of Discrete Dantzig Selector. We mention another MILO formulation for Problem (2), based on
Specially Ordered Sets [7]. We introduce binary variables zi ∈ {0, 1}, which satisfy the condition
(1 − zi )βi = 0 for all i = 1, . . . , p — in other words, if zi = 0, then βi = 0, and if zi = 1, then βi is
unconstrained. This condition can be modeled via integer optimization using Specially Ordered
Sets of Type 1 (SOS-1). More specifically,
(1 − zi )βi = 0

⇐⇒

(βi , 1 − zi ) : SOS-1,

for every i = 1, . . . , p. This leads to the following MILO formulation for Problem (2):
min
β,z

s.t.

p
P

zi

i=1

−δ ≤ dj − hqj , βi ≤ δ

j = 1, . . . , p

(βj , 1 − zj ) : SOS-1

j = 1, . . . , p

zj ∈ {0, 1}

j = 1, . . . , p,

(19)

where, we use the notation as used in Problem (6). Observe that unlike (6), Problem (19)
does not contain any parameter MU in its formulation. Problem (19) may be preferred over
Problem (6) when the different nonzero values of |βbi |’s have widely different amplitudes. In
general, however, we found empirically that the algorithmic performances of formulations (19)
and (6) are comparable. The MILO formulations (6) and (19) are found to work quite well in
obtaining good upper bounds for up to p = 10,000, once they are warm-started via the discrete
first-order methods described in Section 4.1. If additional problem-specific information which
we refer to as “intelligence” is supplied to the MILO formulations (6) and (19), the results are
found to improve substantially — as shown in Section 6.2. More specifically, we use the term
“intelligence” to broadly refer to two components:
(a) providing an advanced warm-start to the MILO solver, obtained via our discrete first-order
methods
(b) arming the MILO solver with information in the form of interval bounds on the regression
coefficients βj , predictions x>
i β, and also bounds on kβk1 and kXβk1 .
We note that the resulting formulation with the additional bounds as suggested in (b), above,
should lead to a solution for Problem (2). We, thus, present the following structured version of
21

formulation (2):
min
β,z

s.t.

p
X

zi

i=1

− δ ≤ dj − hxj , ξi ≤ δ
−

MjU zj

≤ βj ≤

MjU zj

zj ∈ {0, 1}

j = 1, . . . , p

(20a)

j = 1, . . . , p

(20b)

j = 1, . . . , p

ξ = Xβ

(20c)

− MjU ≤ βj ≤ MjU
−

Mξ,i
U

≤ ξi ≤

Mξ,i
U

j = 1, . . . , p

(20d)

i = 1, . . . , n

(20e)

kβk1 ≤ M`
n
X

|ξi | ≤ Mξ` ,

(20f)

i=1
ξ
where the optimization variables are β ∈ Rp , z ∈ {0, 1}p , ξ ∈ Rn , and the parameters MiU , Mξ,i
U , M` , M`
control, respectively, upper bounds on |βi |, |hxi , βi|, kβk1 and kXβk1 . We note that the parameter MU in Problem (6) is such that MU ≥ maxj MjU for j = 1, . . . , p. Problem (20) is equivalent
to the following constrained version of Problem (2):

min

kβk0

s.t.

kX> (y − Xβ)k∞ ≤ δ

β

−MjU ≤ βj ≤ MjU
|hxi , βi| ≤

Mξ,i
U

j = 1, . . . , p

(21)

i = 1, . . . , n

kβk1 ≤ M`
kXβk1 ≤ Mξ` .
Section 5.1 presents several strategies to compute these parameters such that a solution to Problem (21) is also a solution to Problem (2).
We present a few variations of formulation (20) that might be preferred from a computational
viewpoint, depending upon the problem instance under consideration. For large values of p and
n (approximately a few thousand), the constraints appearing in (20a) and (20c) may be replaced
by:
−δ ≤ dj − hx̃j , ξi ≤ δ, ξ = X̃β,
where, X̃> X̃ = X> X and X̃ is triangular — this leads to a sparse representation of the constraints
appearing in (20). When n is large and p is smaller, it may be useful to perform a variable
reduction by removing the variable ξ from (20). This will replace constraint (20a) by −δ ≤
22

dj − hqj , βi ≤ δ, and constraints (20c), (20e) and (20f) will be dropped. Constraints (20b) imply
the bounds indicated in the constraints (20d), hence the constraints (20d) may be dropped in
favor of a formulation with fewer constraints.
We note that formulation (20) is an optimization problem with many more continuous variables
than formulation (6). This implies that the MILO solver needs to do more work at every node,
by solving larger convex linear programs. However, the advantage is that the resulting formulation is more structured, and, thus, tighter lower bounds may be obtained by exploring fewer
nodes. Section 6.2 presents some computational results illustrating the performance of the above
framework.

5.1

Specification of Parameters

We present herein, several data-driven ways to compute the parameters in formulation (20). The
methods presented here are quite different from those proposed in [8], where, the authors rely
crucially on being able to compute analytic expressions for least squares solutions for a given
subset size—such expressions are not available for Problem (2).

5.1.1

Specification of Parameters via Linear Optimization

We present several methods based on linear optimization that can be used to estimate the pab a solution to
rameters appearing in Problem (20), in such a way that these estimates lead to β,
Problem (2).
Bounds on βbi ’s. Consider the following pair of linear optimization problems:
µ+
i :=
s.t.

max βi
β

kX> (y

− Xβ)k∞ ≤ δ,
(22)

µ−
i :=
s.t.

min βi
β

kX> (y

− Xβ)k∞ ≤ δ,

−
+
b
for i = 1, . . . , p. Note that µ+
i and µi provide upper and lower bounds on βi for every i. µi is
typically a strict upper bound to βbi , because (22) does not account for the fact that solutions
b
to Problem (2) are sparse. Similarly, µ−
i is a lower bound to βi , and it is easy to see that
+
−
i
MU = max{|µi |, |µi |} is an upper bound to |βbi |. Note that solutions to Problem (22) are finite
only if the feasible set is bounded. If n > p and if the entries of X are drawn from a continuous
probability measure, then the bounds are finite with probability one. The above bounds can

23

be made tighter by using information about upper bounds on Problem (2) as obtained via the
discrete first-order methods. We describe such methods in Section B.6 (Appendix). Once upper
b ∞ and kβk
b 1
bounds on |βbi |, i.e. MiU , are obtained, they can be used to compute bounds on kβk
as follows:
α0
X
(i)
b ∞ ≤ MU = max Mi and kβk
b 1≤
kβk
MU ,
U
i=1,...,p

i=1
(1)

(2)

(p)

where, α0 denotes an upper bound to Problem (2) and MU ≥ MU ≥ . . . ≥ MU .
b
b can be obtained by solving the following pair of linear
Bounds on hxi , βi’s.
Bounds on hxi , βi
optimization problems:
vi+ (α0 ) :=

max hxi , βi

s.t.

kX> (y − Xβ)k∞ ≤ δ

β

kβk∞ ≤ MU
kβk1 ≤ MU α0
(23)
vi− (α0 ) :=

min hxi , βi

s.t.

kX> (y − Xβ)k∞ ≤ δ

β

kβk∞ ≤ MU
kβk1 ≤ MU α0 ,
for every i = 1, . . . , n.
Analogous to the bounds derived via Problem (22), it is also possible to compute more conservative
b by dropping the constraints kβk∞ ≤ MU and kβk1 ≤ MU α0 in Problem (23).
bounds on hxi , βi,
This gives nontrivial bounds even for the under-determined n < p case, as long as X has rank
n (this is in contrast with the bounds from Problem (22) being vacuous when n < p). It is also
b by including an additional constraint: kXβk∞ ≤ Mξ , and
possible to estimate bounds on hxi , βi
U
using an iterative method as described in the Appendix, Section B.6 (see Step-1–Step-4) while
computing bounds on the regression coefficients.
b
The quantity vi = max{vi+ (α0 ), −vi− (α0 )} provides an upper bound to |hxi , βi|.
In particular,
this leads to the following upper bounds:
b ∞ ≤ max vi
kXβk
i=1,...,n

and

b 1≤
kXβk

n
X

vi ,

i=1

leading to a data-driven method to estimate bounds appearing in (20).
Computational Cost. Computing the quantities appearing in (22), (37) and (23) requires solv-

24

ing at least 2(p + n) linear optimization problems. However, these individual problems are quite
simple to parallelize and they need to be solved once, before proceeding to solve Problem (20).
These linear optimization problems can be solved by simplex based solvers quite easily for p in
the lower thousands (typically less than a minute with Gurobi’s simplex solver).
5.1.2

Specification of Parameters from warm-starts

We present herein, simple practical methods to compute the parameter values by using good
b 0 denote a solution that corresponds to a good upper bound
upper bounds to Problem (2). Let β
b 0 can be obtained from Algorithm 3, for example. One can also use the
to Problem (2). β
solution obtained from Algorithm 3 as a warm-start to Problem (19) and allow it to run for
b 0 . The parameters appearing in the
a few minutes — the resulting estimate may be used as β
b 0 , as follows. To be on the conservative side, we recommend setting
bounds can be based on β
b 0 k∞ .
the same bound for all the MiU ’s: for example, they can all be assigned the value τ kβ
b0
Similarly, anconservative bound foroall the Mξ,i
U ’s is given by τ kXβ k∞ . In addition, we can set
b 0 k∞ , for some value of τ ∈ {1.5, 2}.
b 0 k0 kβ
b 0 k∞ , τ kβ
b 0 k1 and Mξ = τ kXβ
M` = min τ kβ
`
The method described above leads to parameter specific bounds as a simple by-product of our
general algorithmic framework. Unlike the methods in Section 5.1.1, it requires no additional
computation. On the other hand, the bounds in Section 5.1.1 are conservative, because they are
implied by the bounds from Problem (2).

6

Numerical Experiments: Algorithmic Performance

In this section, we report extensive numerical experiments that demonstrate: (a) the usefulness of
the discrete first-order methods (Section 4.1) in obtaining good quality upper bounds, especially
when they are used to provide warm-starts to MILO solvers — this is shown in Section 6.1; and
(b) how advanced warm-starts, coupled with the enhanced formulations presented in Section 5.1,
can be used to improve the overall run-time for off-the-shelf MILO solvers, when proving global
optimality for the Discrete Dantzig Selector problem — this is shown in Section 6.2.
All computations were carried out on Columbia University’s high performance computing (HPC)
facility, http://hpc.cc.columbia.edu/, on the Yeti cluster computing environment. The discrete first-order methods were implemented in Matlab 2014a, and we used Gurobi [25] version
6.0.3. For all experiments in Sections 6.1 and 6.2 (except the large scale examples) we used 16GB
of memory.

25

Type-1 (n = 100, p = 1000)
Data Fidelity
Parameter
1.5δ̄

δ̄

0.5δ̄

0.2δ̄

Time
(in secs)
2 (*)
120
500
90 (*)
120
500
57 (*)
120
500
120
210 (*)
500

Type-2 (n = 300, p = 1000)

Quality of Upper Bounds
With Warm
Vanilla
0
–
0
0
0
0
0
42.85
0
42.85
0
28.57
0
200
0
86.66
0
26.66
3.12
25
0
15.62
0
15.62

Data Fidelity
Parameter
1.5δ̄

δ̄

0.5δ̄

0.2δ̄

Type-3 (n = 600, p = 2000)
Data Fidelity
Parameter
1.5δ̄

δ̄

0.5δ̄

0.2δ̄

Time
(in secs)
500
530 (*)
950
500
875 (*)
950
55 (*)
500
950
500
560 (*)
950

Time
(in secs)
50 (*)
120
500
120
132(*)
500
35 (*)
120
500
40 (*)
120
500

Quality of Upper Bounds
With Warm
Vanilla
0
–
0
214.28
0
14.28
6.66
146.66
0
73.33
0
0
0
–
0
29.62
0
25.92
0
–
0
73.44
0
23.43

Type-4 (n = 58, p = 2000)

Quality of Upper Bounds
With Warm
Vanilla
7.14
–
0
–
0
28.57
11.11
–
0
137.03
0
33.33
0
–
0
–
0
120.37
0.8
–
0
–
0
77.6

Data Fidelity
Parameter
1.2δ̄

δ̄

0.5δ̄

0.2δ̄

Time
(in secs)
300 (*)
370
600
300
367 (*)
600
300
560 (*)
600
145 (*)
300
600

Quality of Upper Bounds
With Warm
Vanilla
0
220
0
0
0
0
16.66
–
0
216.66
0
0
5
–
0
95
0
95
0
–
0
165
0
60

Table 1: Tables showing “Quality of Upper Bounds”, defined as 100 × (halg − b
h)/b
h, where halg refers to
b
the objective value obtained by algorithm “alg” (at the given time), and h is the best objective value found
in the entire run-time duration of all the algorithms. Two cases of “alg” ∈ {“With Warm”, “Vanilla”}
have been considered: “With Warm” denotes MILO warm-started with a solution from Algorithm 3 (the
structured formulations of Section 5 are not used here); “Vanilla” denotes a MILO solver without any
warm-start specification. “With Warm” is found to obtain the best upper bound for a given computation
time-limit in all the instances. In several instances, MILO is found to improve the solution obtained via
Algorithm 3, after accepting it as a warm-start. For Type-1,2 the total time limit was 500 secs; for Type-3
it was 950 secs, and for Type-4 the algorithms were considered for a total time limit of 600 secs. For method
“With Warm”, the times reported show the overall time taken by Algorithm 3 and the MILO algorithm.
An asterisk “(*)” indicates that the best solution is obtained at that time. A “–” means that no feasible
solution was obtained by the algorithm in that time.

26

6.1

Obtaining Good Quality Upper Bounds

From a practical viewpoint, being able to obtain good quality upper bounds to Problem (2) is,
perhaps, of foremost importance. To demonstrate the effectiveness of our computational framework in this regard, we perform a series of experiments on the data-types described below.
Type-Synth: We generate a Gaussian ensemble Xn×p ∼ MVN(0, Σ), where σij = ρ|i−j| for some
value of ρ ∈ [0, 1), with the convention that 00 = 1. The underlying true regression coefficient
vector, β ∗ ∈ Rp , has βj∗ = 1 for k ∗ equi-spaced values of j ∈ {1, . . . , p} and βj∗ = 0 for the
remaining values of j.
Type-1: This is of Type-Synth with n = 100, p = 1000, ρ = 0, k ∗ = 10. We studied
Problem (2) for four different values of the parameter δ set at δ̄(1.5, 1, 0.5, 0.2) with δ̄ being
defined below.
Type-2: This is of Type-Synth with n = 300, p = 1000, ρ = 0.8, k ∗ = 25. Here δ values were
set as δ̄(1.5, 1, 0.5, 0.2).
Type-3: This is of Type-Synth with n = 600, p = 2000, ρ = 0.8, k ∗ = 40. Here δ values were
set as δ̄(1.5, 1, 0.5, 0.2).
Type-4: This is a semi-synthetic dataset: we
pression dataset7 from Ch. 16 of the book [26].
p = 2000 and there were n = 58 observations.
model with kβ ∗ k0 = 10, βj∗ = 1 for j ≤ 10, and
δ̄(1.2, 1, 0.5, 0.2).

considered the Radiation sensitivity gene exThe features were randomly downsampled to
We generated response y based on a linear
βj0 = 0 for j > 10. Here δ values were set as

Type-5: This is of Type-Synth with n = 1000, p = 3000, ρ = 0, k ∗ = 10; we considered one
value of δ, which was set to δ.
In each of the above examples, after X was generated, we standardized its columns to have zero
iid
mean and unit `2 -norm. Then, the response was generated as y = Xβ ∗ + , where i ∼ N (0, σ 2 ),
and σ 2 was adjusted to match the selected value of SNR= Var(x> β ∗ )/σ 2 (taken as 3 in all the
above cases); the reference value of the tuning parameter was set to δ̄ = kX> (y − Xβ ∗ )k∞ .
We studied different first-order algorithms described in Section 4. Algorithm 3 was empirically
seen to have the best performance over its constituents, Algorithms 1 and 2, when used separately.
Hence, we used Algorithm 3 in all the experiments to obtain good upper bounds to Problem (2).
The solution obtained from Algorithm 3 was passed as a warm-start to the MILO formulation (6)
(for a large value of MU = 103 ) — this hybrid MILO approach is denoted by “With Warm” in
Table 1. We compared this method with the vanilla MILO formulation (6) (for a large value of
MU = 103 ), which was implemented without any warm-start information. Table 1 shows the
7

We downloaded the dataset from the website http://statweb.stanford.edu/~tibs/ElemStatLearn/
datasets/

27

objective values obtained by these two methods — the MILO algorithm aided with advanced warmstarts was found to perform the best across all the examples. For the hybrid approach (“With
Warm”), in many of the instances, the solution obtained by Algorithm 3 was further improved by
MILO. In some cases, the vanilla MILO approach took a while before it was able to find a feasible
solution. For example, in the Type-5 setting (which does not appear in Table 1) the best solution
was delivered by Algorithm 3 within one minute; in contrast, the vanilla MILO algorithm failed to
find a feasible solution within 1000 seconds.

●
●
●
●
●
●
●

●
●

●
●

Vanilla
With Intelligence

●

●
●

60
●
●
●
●

40

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

500

●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●

20

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●●
●

●

●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●

20

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

60

●

40

60
40

●

●
●

●
●

●

●
●
●

●
●
●
●
●
●
●

0

20

Vanilla
With Intelligence

●

●

0

●

●
●
●
●
●

80

●

●

100

Vanilla
With Intelligence

●
●
●●
●
●
●
●

0

MILO Optimality Gap (in %)

80

●

80

●

b 0 =7
kβk

1000

1500

2000

Time (secs)

2500

3000

0

100

●

100

Diabetes Dataset (n = 442, p = 64)
b 0 = 14
kβk

b 0 = 31
kβk

● ●

0

500

1000

1500

2000

2500

3000

Time (secs)

●

0

●

500

1000

1500

2000

2500

3000

Time (secs)

Figure 3: The evolution of MILO Optimality gaps (defined in Figure 2) as functions of time (in secs) for
the MILO methods with and without problem-specific information. We consider three different values of
the data fidelity parameter δ, leading to solutions with different number of nonzeros as mentioned in the
figure panels. Here, “With Intelligence” refers to a MILO algorithm aided with an advanced warm-start, in
addition to the bounds specified in Section 5, and “Vanilla” refers to a MILO algorithm without any such
additional information. MILO is found to benefit significantly from additional problem-specific information.

6.2

Lower Bounds and Certificates of Optimality

Here we demonstrate how our framework delivers certifiably optimal solutions to Problem (2).
In our first series of experiments we considered the popular diabetes dataset [19], which we
examined with interaction terms included, giving us p = 64 and n = 442. All the features and
the response were mean-centered and standardized to have unit `2 -norm. Figure 3 shows the
performance of two versions of MILO – “With Intelligence” and “Vanilla”. “With Intelligence”
refers to MILO formulation (20), where a MILO solver is provided with an advanced warm-start,
b 0 . The parameter specifications are obtained based on the method in Section 5.1.2; here, we
say, β
used the box constraints and `1 -constraint on β. The “Vanilla” version of MILO was not provided
with any such problem-specific information — we used formulation (6), as in Section 6.1. Our
experimental results (Figure 3) show that “Intelligence” significantly enhances the performance
of the MILO solver, in terms of proving global optimality. Usually, we observe that for a fixed n, p
with n > p the time to certify optimality is smaller when k is small or close to p – intuitively, this
28

100

With Intelligence

●

With Intelligence

80

80
60
●

●● ● ●
●

●
● ● ●

●
●● ● ● ● ● ●

20

20

●●

20

● ● ●

●

●

●

●

●

●

● ● ● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●
●● ●

20

● ● ● ● ●●

40

●●● ●

●

b 0 ) = (900, 2500, 20)
(n, p, kβk

60

80

100

100

b 0 ) = (750, 2000, 20)
(n, p, kβk

40

●●●● ●●

With Intelligence

●

40

60

●

b 0 ) = (600, 1500, 15)
(n, p, kβk

60

100

With Intelligence

80

●

40

MILO Optimality Gap (in %)

b 0 ) = (400, 1000, 20)
(n, p, kβk

● ●
●

●● ●

●

●●
●●●● ● ● ●

●

● ●
● ●

●● ●

200

400

600

800

Time (secs)

1000

100

150

200

● ●

250

●

200

Time (secs)

300

400

Time (secs)

500

0

●
●●

0

●

●

0

0

●

●

600

800

1000 1200 1400 1600

Time (secs)

Figure 4: The evolution of MILO Optimality gaps (defined in Figure 2) as a function of time (after a warmstart was supplied to it) for different synthetic examples with varying problem-sizes. “With Intelligence”
refers to MILO aided with an advanced warm-start information, in addition to the bounds specified in
Section 5, as described in the text. The “Vanilla” version of the MILO was found to perform worse, and is,
thus, not shown in the figures. In all these instances, the `1 -Dantzig Selector resulted in denser solutions.

is due to the “search-space” being small. The computation time increases as k becomes closer to
p/2. This is reflected in Figures 2 and 3.

6.2.1

Moderate Scale Examples

We considered some examples of Type-Synth for ρ = 0 and different values of n, p, k ∗ ; in all the
examples, we set δ = kX> (Y − Xβ ∗ )k∞ . The results for MILO with intelligence are displayed
b 0 ) from a combination of Algorithm 3 and
in Figure 4. We obtained an advanced warm-start (β
MILO formulation (19), where the latter was allowed to run for an overall time limit of 500 seconds.
b 0 was used to initialize formulation (20) — the parameter specifications in the
The warm start β
formulation were obtained based on the method in Section 5.1.2. We also experimented with the
version of formulation (20) that considers only box constraints on β; the results were often found
to be roughly similar — both methods certified optimality, though there were some differences
in the total run-time (roughly around a few minutes). For all the synthetic examples presented
in Figure 4, the vanilla version of MILO took much longer to prove optimality and hence they are
not shown in Figure 4.
In all these instances, the `1 -Dantzig Selector, not surprisingly, resulted in a solution that was
more dense than the corresponding Discrete Dantzig Selector.

29

n
2,000
4,000
7,000
6,000
4,000
3,000
3,000
4,000
1,000
5,000
5,000
10,000

p
5,000
5,000
7,000
7,000
8,000
8,000
9,000
9,000
10,000
10,000
10,000
10,000

k∗
30
60
20
60
20
20
20
20
10
10
20
30

(Synthetic Examples)
Upper Bound Lower Bound
30
30
60
60
20
20
60
60
20
20
20
20
20
20
20
20
10
10
10
10
20
20
30
27

MILO Gap
0
0
0
0
0
0
0
0
0
0
0
10%

Time (hrs)
8.3
20.0
21.7
20.0
41.9
18.3
47.2
44.4
14.2
2.5
47.5
42.5

n
6,000
6,000

p
4,500
4,500

k∗
20
40

(Real Data Examples)
Upper Bound Lower Bound
20
20
40
37

MILO Gap
0
10%

Time (hrs)
5.0
12.5

Table 2: Table showing times taken to reach global optimal solutions for several problem instances (both
synthetic and real data), with p up to 10,000. For each instance, we list the upper bound, the lower
bound, the corresponding MILO (Optimality) Gap and the time taken to reach the listed lower bound by
the MILO solver equipped with “Intelligence” (as described in Figure 4). The times (in hours) refer to
those taken by the MILO solver after being provided with a warm-start. In all the instances, the best
upper bound was obtained around approximately one hour, however, it took much longer to obtain a
certificate of global optimality via the (almost) matching lower bounds. The MILO Gap was found to be
zero in all the instances apart from the two where the algorithm was terminated upon obtaining a lower
bound within 10% of the upper bound. The results demonstrate that certifiably optimal Discrete Dantzig
Selector solutions can be obtained for large scale instances.

6.2.2

Large Scale Examples

We consider several large scale examples with p ranging in 4,500 to 10,000: these problem-sizes
are orders of magnitude greater than those considered in [8]. These computations were performed
with 100 GB memory.
We studied a host of synthetic examples, all generated as in Section 6.2.1. We also considered a
semi-synthetic dataset derived from the well-known Gisette data http://archive.ics.uci.edu/
ml/datasets/Gisette. Here, we generated a response y, based on the Gisette data covariates
(each feature was standardized to have zero mean and unit `2 norm) — here, n = 6000 and
p = 4500, we set βi∗ = 1, i = 1, . . . , k ∗ and the remaining βi∗ = 0; SNR=3 and considered two
instances with k ∗ = 20, 30.
The algorithmic set-up was similar to that used in Section 6.2.1. The results are presented in
Table 2. In all these instances, the `1 -Dantzig Selector resulted in a solution that was more
dense than those obtained via the Discrete Dantzig Selector. For all the synthetic examples,
the MILO solver delivered solutions that matched the optimal solution of the data generating
mechanism. Typically, the time taken to prove optimality marginally increases with larger values
of k ∗ (for a fixed n, p); for a fixed k ∗ , p the times taken to certify optimality increases with
decreasing values of n. The examples demonstrated in this paper, show the largest instances
30

of discrete optimization problems for exact variable selection, that can be solved to provable
optimality.

7

Numerical Experiments: Statistical Properties

We conducted a series of synthetic experiments to understand the statistical properties of the Discrete Dantzig Selector and compare them to those of the `1 -Dantzig Selector and variants.
We used the following datasets in our analysis.
Example-A: This is of Type-Synth with n = 200, p = 500, ρ = 0 and k ∗ = 20.
Example-B: This dataset was similar to the one taken in Example-A, but the amplitudes and
signs of the true regression coefficients were allowed to vary: the twenty nonzero βj∗ ’s were equally
spaced in the interval [−10, 10].
Example-C: This is of Type-Synth with n = 100, p = 500, ρ = 0.85, k ∗ = 10.
Example-D: We set n = 100, p = 300 and let X ∼ MVN(0, Σ), where σ12 = σ21 = 0.7, and
all the remaining σjk are equal to zero. We also took β1∗ = 1 and β2∗ = −1, with the remaining
coefficients βj∗ set to zero, resulting in k ∗ = 2. (This example is a larger version of Example 10
described in Section 1 and illustrated in Figure 1.)
In each of the above cases, after X was generated, we standardized its columns to have unit
iid
`2 -norm. Then, the response was generated as y = Xβ ∗ + , where i ∼ N (0, σ 2 ), and σ 2 was adjusted to match the selected value of SNR, which was varied across {3, 10} in the examples.
We considered the following estimators in our analysis:
• “Warm” — this method applies a heuristic strategy to obtain upper bounds to Problem (2).
We used8 Algorithm 2, described in Section 4.1.2.
• “L0-DS” — the solution obtained from “Warm” is taken as a warm-start to a MILO solver
and subsequently allowed to run with a time limit of 4000 seconds.
• “L0-DS-Pol” — this is a “polished” version of the Discrete Dantzig Selector estimator “L0DS” and is obtained by performing a simple least squares fit on the support of the “L0-DS”
estimate.
• “L1-DS” — this is the original `1 -Dantzig Selector.
• “L1-DS-Pol” — this is a polished version of the “L1-DS”.
8

This is similar to a re-weighted `1 -minimization [15] method applied to Problem (2). We took the penalty
ργ (|β|) ∝ log(|β|/γ + 1) on a geometrically decreasing grid of ten γ values: γi = 10−2 × 0.8i−1 for i = 1, . . . , 10.

31

Each of the above estimators were computed on a range of approximately thirty different δ values
around δ̄ = kX> (y − X> β ∗ )k∞ . We considered ten different replications (based on different 
realizations) and took the median of the results. The optimal tuning parameter (δ opt ) for every
model was selected based on the value of δ that minimized the estimation error with respect to the
true regression coefficients. For this chosen value of δ opt we considered different metrics to assess
the performance of the different estimators. We computed the squared `2 -error in estimating the
b − β ∗ k2 . We also considered the “Variable Selection error”, which is
regression coefficients: kβ
2
Pp
∗
b and S ∗ is the jth
b := Supp(β),
b
defined as j=1 1(Sj 6= Sj ), where Sbj is the jth coordinate of S
j
element of S∗ = Supp(β ∗ ). Finally, we computed the “number of nonzeros”, which refers to the
b
number of nonzero coefficients in β.
A collection of representative results with SNR=10 is displayed in Figure 5. The error bars
√
correspond to standard errors, the width being set to 2ŝ/ N where, ŝ is the mean absolute
deviation around the median and N denotes the number of replicates (here, ten). A larger display
of additional examples with varying SNR values is presented in Table 4 in Appendix C, where we
b − Xβ ∗ k2 /kXβ ∗ k2 . In Table 5, given in the
also report the “Prediction Error”, defined as kXβ
2
2
same section, we provide comparisons with the polished version of the `1 -Dantzig Selector. Our
experiments show that polishing the `1 -Dantzig Selector may lead to marginally better solutions
relative to the original Dantzig Selector, but the corresponding statistical performance is inferior
to that of the estimates based on the Discrete Dantzig Selector estimator.
We note that the performance of the Lasso was found to be quite similar to that of the `1 -Dantzig
Selector. The statistical performance of the subset selection procedure (4), as described in [8], was
found to be similar to that of the Discrete Dantzig Selector for p ≤ 1000. Because the main focus
of the paper is to show that Discrete Dantzig Selector is a computationally tractable procedure,
which delivers estimates with better statistical properties than its `1 counterpart, we restrict our
numerical studies to the methods listed above.
Summary of Findings. Based on the experimental results, we observe that the Discrete Dantzig
Selector and its polished variant perform quite well, when compared to the competing methods
in terms of estimating β ∗ ; they also demonstrate superior variable selection properties (not surprisingly, the `0 methods obtain the sparsest models across all the examples). “Warm” does not
perform very well when compared to “L0-DS”, even though both methods attempt to solve Problem (2) — this suggests that estimators based on rigorous optimization procedures have better
statistical properties. We observe that “L0-DS” and “L0-DS-Pol” possess similar variable selection properties, however, the latter may lead to better estimators of β ∗ and Xβ ∗ , due to the least
squares post-processing. Polishing of the `1 -Dantzig Selector may not lead to better solutions,
due to the weak variable selection properties of “L1-DS”. In some cases, when the value of ρ is
quite large and, consequently, the covariates are highly correlated (see Example-C, Example-D),
the basic problem of variable selection becomes difficult: instead of choosing a “signal” variable,
the “L0-DS” chooses its correlated surrogate. In these cases, as expected, we observe that the
32

Example-A

Example-C

Example-D

(n =Indn_p_snr_2005003−snr3−bet
200, p = 500)

(n =Decn_p_snr_1005003−snr3−bet
100, p = 500)

(n =
100, p = 300)
Petn_p_snr_1003003−snr3−bet

b − β ∗ k2
kβ
2

1.5

0.06

3.0

0.05

2.5
1.0

0.5

0.0

2.0

0.04

1.5

0.03

1.0

0.02

0.5

0.01

0.0
L0−DS L0−DS−Pol Warm

L1−DS

0.00
L0−DS L0−DS−Pol Warm

Indn_p_snr_2005003−snr3−err01

L1−DS

L0−DS L0−DS−Pol Warm

Decn_p_snr_1005003−snr3−err01

L1−DS

Petn_p_snr_1003003−snr3−err01

40

Error in Variable Selection

100

35
30

30

80

25
60

20

20

15

40

10

10
20

5
0

0
L0−DS L0−DS−Pol Warm

L1−DS

0
L0−DS L0−DS−Pol Warm

Indn_p_snr_2005003−snr3−nnz

L1−DS

L0−DS L0−DS−Pol Warm

Decn_p_snr_1005003−snr3−nnz

L1−DS

Petn_p_snr_1003003−snr3−nnz

40

50

Number of nonzeros

120
40

100
80

30

30
20

60
20
40

10
10

20
0

0
L0−DS L0−DS−Pol Warm

L1−DS

0
L0−DS L0−DS−Pol Warm

L1−DS

L0−DS L0−DS−Pol Warm

L1−DS

Figure 5: The statistical performance of the Discrete Dantzig Selector (L0-DS); its polished version,
which does a least squares refitting on the support (L0-DS-Pol); the heuristic estimates delivered by
Algorithm 2 (Warm); and the original `1 -Dantzig Selector (L1-DS). We display three different metrics:
[top panel] squared `2 -error in estimating β, [middle panel] the 0-1 variable selection error and [bottom
panel] the number of nonzeros in the optimally selected model; the horizontal dotted line shows the number
of nonzeros in the “true” model. We observe that the Discrete Dantzig Selector based approaches perform
very well in terms of obtaining a model with high quality estimation and variable selection properties. Their
models are substantially sparser than those for the `1 -based methods. The heuristic approach, “Warm”
which approximately optimizes Problem (2), falls short in obtaining high quality statistical estimates. A
number of additional experiments (with results similar to this figure), are presented in Section C of the
Appendix.

33

Example
ρ=0
n = 100
p = 200
k ∗ = 20
ρ = 0.7
n = 100
p = 200
k ∗ = 10
ρ = 0.7
n = 300
p = 600
k ∗ = 20
ρ = 0.7
n = 300
p = 1000
k ∗ = 25

Metric
Variable Selection
Error
b − β ∗ k2
kβ
2
Variable Selection
Error
b − β ∗ k2
kβ
2
Variable Selection
Error
b − β ∗ k2
kβ
2
Variable Selection
Error
b − β ∗ k2
kβ
2

Time (secs)
60
200
60
200
60
200
60
200
100
500
100
500
150
1000
150
1000

L0-DS-Pol
0
0
0.309
0.309
0
0
0.063
0.063
0
0
0.094
0.094
0
0
0.175
0.175

MIPGO
20
14
20.00
8.79
10
4
0.271
0.076
9
5
9.142
1.548
25
16
25.00
14.893

Table 3: In all the above instances, the generated data is of Type-Synth, with SNR=10. Both the
Discrete Dantzig Selector and MIPGO (with the MCP penalty) methods were run on twenty different
values of the tuning parameter, and the best solutions are reported. L0-DS-Pol refers to the least squares
solution obtained on the variables selected by Discrete Dantzig Selector. For every value of the tuning
parameter, each method was run for a time budget of t1 , t2 seconds with t1 < t2 specified in the “Time”
column. The Discrete Dantzig Selector method reaches the best solution within t1 seconds, much earlier
than its competitor. The MIPGO method is seen to take orders of magnitude longer to get a solution of the
same quality as the Discrete Dantzig Selector — the differences become more pronounced with increasing
problem size.

“L0-DS” incurs relatively large variable selection error–the prediction accuracy of these models however, demonstrate a more optimistic picture than the variable selection properties (See
Table 4).

7.1

Comparisons with Least Squares Subset Selection

We discuss some comparisons of our proposal with the recently proposed methods: Problem (4)
by [8] and MIPGO [33].
For underdetermined problems (with n < p) the authors in [8] (see Section 5.3.2 in [8]) point
out that MIQO solvers take a long time to certify optimality, by producing matching upper and
lower bounds. For problems with n ≤ 100 and p = 1000, [8] demonstrated how the MIQO methods
for Problem (4) could certify local optimality9 in a (small) bounding box around a candidate
solution. We observed in our computational experiments that Problem (2) is orders of magnitude
faster than (4) for underdetermined problems, in obtaining solutions with certificates of global
9

We note that certifying local optimality i.e., optimality in a neighborhood of a candidate solution is also an
NP-hard problem.

34

optimality. On several randomly generated problem instances generated as per Type-Synth with
n = 100, p = 200, ρ = 0, k ∗ = 10 and SNR=10, Problem (2) was solved to global optimality,
i.e., zero optimality gap with a median time of about 4 minutes. On the same instances, the
MIQO formulation for Problem (4) took more than 7 hours of computation time to obtain similar
optimality certificates. In addition, a MILO formulation for Problem (2) consumes much less
memory than a comparable MIQO formulation for Problem (4). For example, on problem instances
with n = p ∈ {1.5, 2, 2.5, 3} × 103 , we observed that Problem (4) requires at least twice as much
memory as that for Problem (2), within the first 800 seconds of computation time. The memory
requirement for a MIQO for Problem (4) with n = p = 3000 was more than 12GB.
MIPGO [33] is a discrete optimization framework for minimizing a regularized version of the least
squares loss, with a nonconvex quadratic penalty (for example, SCAD or MCP). This corresponds
to a nonconvex quadratic optimization problem, which the authors express as a discrete linear
optimization problem via linear complementary constraints [24]. This representation results in
many more binary variables (several multiples of p) than that required for the Discrete Dantzig
Selector. For example, in the case of the MCP penalty, the paper [33] presents a MILO with 4p
binary variables and many more continuous variables. In particular, with n = 300 and p = 1000,
the MIPGO solver10 creates a problem with 38,000 variables and 27,000 equality constraints,
which after presolve reduces to a problem with approximately 7,000 continuous and 4,000 binary
variables. These optimization problems are substantially larger than (6), which is a MILO with p
continuous and as many binary variables. As a result, the MIPGO formulation seems to become
computationally expensive as the dimensionality of the problem increases. This is illustrated in
Table 3 – we show on several synthetic instances that, with a particular time budget, the Discrete
Dantzig Selector formulation (6) equipped with a warm-start from Algorithm 2, leads to better
solutions than MIPGO. The quality of solutions produced by MIPGO is found to improve with
more computation time, but the time taken can be substantially larger than that of the Discrete
Dantzig Selector. The synthetic datasets for Table 3 were generated the same way as in Section 6.
We set the concavity parameter for the MCP penalty in the MIPGO code to its default value of
a = 2.

Acknowledgements
The authors thank the anonymous referees for their helpful comments that led to improvements in
the manuscript. The authors also thank Emmanuel Candes and Robert Freund for helpful suggestions and encouragement. R.M. thanks Jonathan Goetz, Juan-Pablo Vielma and Dimitris Bertsimas for helpful discussions. R.M.’s research was partially supported by ONR N000141512342 and
a grant from the Moore Sloan Foundation. Peter Radchenko’s research was partially supported
10
We used the code of [33], obtained from the first author’s website. The numbers are read off from the Gurobi
log report.

35

by NSF Grant DMS-1209057.

A

Proofs for Section 3

Proof of Theorem 1. Let ˜ denote the projection of y onto the orthogonal complement to the
space spanned by the predictor vectors xj , j = 1, ...p. Note that
ky − Xβk2 =

p
X
(cj − βj )2 + k˜
k2 .
j=1

Thus, under the constraint kβk0 ≤ k, the smallest sum of squares is achieved by setting βj = cj
for j = 1, ..., k and βj = 0 for j = k + 1, ..., p. This completes the proof of part 1.
Note that
|x>
j (y − Xβ)| = |cj − βj |.
Thus, the constraint kX> (y − Xβ)k∞ ≤ δ is satisfied if and only if βj ∈ [cj − δ, cj + δ] for
j = 1, ..., p. In order to minimize kβk0 , coefficients βj for which 0 ∈ [cj − δ, cj + δ] are set to zero.
Thus, βj = 0 if and only if |cj | ≤ δ, which implies βj = 0 for j > k. This completes the proof of
part 2.

Proof of Theorem 2. Throughout this proof we omit the words “with probability tending to
one” to improve the presentation. The constraint kX> (y − Xβ)k∞ ≤ δ implies kX> X(β ∗ − β) +
X> k∞ ≤ δ. Recall that δ = o(n1/2 ) and note that kX> k∞ = Op (1), due to the scaling of the
predictors and the assumptions on the . Consequently, kX> X(β ∗ − β)k2 = op (n1/2 ). Because
X> X converges to an invertible matrix C, we conclude that there exists a op (1) sequence of
random variables bn , such that the bound
b − n−1/2 β ∗ k2 ≤ bn
kn−1/2 β
b Recall that β ∗ = n1/2 β̃ ∗ , for some
simultaneously holds for all the Discrete Dantzig solutions β.
∗
fixed vector β̃ . Consequently, βj∗ 6= 0 implies βbj 6= 0 for j = 1, ..., p. In other words the support
of each Discrete Dantzig solution contains J ∗ . It is only left to show that the cardinality of
each such support cannot be greater than |J ∗ |. Note that with probability tending to one, β ∗ is
feasible for the Discrete Dantzig optimization problem. Indeed,
kX> (y − Xβ ∗ )k∞ = kX> k∞ = Op (1),
b 0 ≤ |J ∗ | holds
which is bounded above by δ, due to the assumption δ → ∞. Thus, inequality kβk
for each Dantzig Selector solution, which completes the proof of part 1.

36

In the paragraph above we deduced that the minimum value of the Discrete Dantzig objective
function equals |J ∗ |. We also showed that β ∗ is feasible for the Discrete Dantzig optimization
problem. A similar argument establishes the feasibility of β O . Consequently, β ∗ and β O are
indeed Discrete Dantzig solutions, which completes the proof of part 2.

Proof of Proposition 1. Consider an arbitrary nonzero θ ∈ Rp , such that kθk0 ≤ 2k. Let J0
be the index set of the k largest, in magnitude, coordinates of θ. Observe that kθ J0c k1 ≤ kθ J0 k1
and kθk2 = kθ J01 k2 , because m ≥ k. Thus
kXθk2
kXθk2
≥ κ(k, c0 , m),
=
kθk2
kθ J01 k2
for c0 ≥ 1, which implies γ(2k) ≥ κ(k, c0 , m). Also note that
2

kXθk22
kXθk22
≥
≥ [κ(k, c0 )]2 ,
kθk22
kθ J0 k22

√
for c0 ≥ 1, which gives γ(2k) ≥ κ(k, c0 )/ 2.
Proof of Theorem 3 and Corollary 2. Note that X>  is a mean zero Gaussian vector, such
that the variance of each component is σ 2 . Consequently, it follows from well-known maximal
inequalities for Gaussian variables that the bound kX> k∞ ≤ δ holds with probability at least
√
1 − (pa π log p)−1 . The rest of the proof is conducted on the set where the above bound is valid.
Note that on this set β ∗ is a feasible solution for the optimization problem (2), which implies
b 0 ≤ kβ ∗ k0 . Recall that we denote kβ ∗ k0 by s∗ and derive the following inequalities:
kβk
b − β∗
γ(2s∗ )2 β

2
2

b − β∗ )
≤ X(β

2
2

b − β ∗ )> X> X(β
b − β∗ )
= (β
b − β∗ )
X> X(β

b
≤
X> (y − Xβ)
≤

∞
∞

b − β∗
β

1

+ n−1 X> 


∞

b − β∗
β

1

.

b ∞ and kX> k∞ are bounded above by δ, we derive
Because both kX> (y − Xβ)k
b − β∗
β

2
2

b − β∗
≤ γ(2s∗ )−2 2δ β

1

.

b − β ∗ k1 ≤ (2s∗ )1/2 kβ
b − β ∗ k2 to either the left or the right hand side of
Applying inequality kβ
the above display yields the `1 and the `2 estimation bounds, respectively, in the statement of
Theorem 3.

37

Finally, to establish the prediction error bound, observe the following inequality:
b − β∗ )
X(β

2
2

b − β∗
≤ 2δ β

1

,

which is is a direct consequence of the two displays given above. We then complete the proof of
Theorem 3 by combining the above display with the inequalities
b − β∗
β

b − β∗
≤ (2s∗ )1/2 β

1

2

b − β∗ )
≤ (2s∗ )1/2 γ(2s∗ )−1 X(β

2

.

Corollary 2 follows directly from the `1 estimation bound in Theorem 3.
Proof of Theorem 4. We again focus on the set of high probability, where inequality kX> k∞ ≤
δ holds. Because β ∗ is a feasible solution to the optimization problem (2), we have kβ ∗ k0 ≥ sbLB .
b 0 ≤ (1 + ψ)kβ ∗ k0 . The rest of the proof is identical to the one for Theorem 3, with one
Thus, kβk
b 0 ≤ kβ ∗ k0 contains an additional factor (1 + ψ).
exception: the bound kβk

B

Additional Algorithm Details and Proofs

B.1

Details on Algorithm 1

Update (11) in the ADMM algorithm can be performed via a hard thresholding operation [18],
as it is of the form:
b 0 ) := arg min kβ − ck2 + λ0 kβk0 ,
β(λ
(24)
2
β

for an appropriately chosen λ0 = λ2 and c = α + λ1 ν; a solution is given by βbj (λ0 ) = cj 1(|cj | >
√
λ0 ), j = 1, . . . , p. The update step (12) involves the following projection:
min f (α) := kα − ck22
α

s.t. kX> (y − Xα)k∞ ≤ δ,

(25)

where c = β − ν/λ. While the projection (25) can be computed by using standard quadratic
programming methods, in our experience, we found them11 to be quite time consuming for larger
problems (p ≥ 1000), especially because this projection needs to be computed for every iteration
(indexed by k) of (11)–(13). Thus, we recommend using specialized first-order methods — these
methods also naturally make use of warm-start information, which is particularly useful to us
due to the iterative nature of the updates (11)–(13). Unless X has uncorrelated columns, it is
not straightforward to solve (25) in its primal form — we thus consider a dual of Problem (25),
11

Our reference is Gurobi’s quadratic programming solver.

38

for which we apply first-order methods for convex composite minimization [36]. To improve the
flow of presentation, we relegate the description of a more general first-order method, which also
applies to Problem (25), to Section B.3 in the Appendix. We repeat steps (11)–(13) until an (approximate) convergence criterion is met — see for example [34, 12] for convergence results for the
general method. We terminate the algorithm as soon as the successive changes in the β updates
become small and one has approximate primal feasibility (see Step (3) in Algorithm 1).

B.2

Additional Details on Algorithm 2: Solving Problem (17)

Observe that Problem (17) is of the composite form [36]:
min f1 (θ) + f2 (θ) s.t. θ ∈ C,
θ

(26)

where the function f1 (θ) is smooth, with its gradient Lipschitz continuous: k∇f1 (θ)−∇f1 (θ 0 )k ≤
Lkθ − θ 0 k; f2 (θ) is nonsmooth and C is a convex set. In our specific case, the smooth component
P
is the zero function, f2 (θ) = pi=1 wi |θi | and C = {θ : kX> (y − Xθ)k∞ ≤ δ}. Thus, one may
appeal to first-order optimization methods [2, 37, 36] for composite function minimization. This
requires solving, at every iteration, a problem of the form:
θ m+1 ∈ arg min
θ

s.t.

L
2 kθ

kX> (y

m

m

− θ k22 +

p
P

wi |θi |

i=1

(27)

− Xθ)k∞ ≤ δ,
m

for some choice of L > 0 and θ , and wi = |ρ0γ (|βik |)|. If θ = θ m , then the above update sequence
becomes identical to proximal gradient descent [2]. One may also use accelerated gradient descent
methods, with a momentum term. We describe in Section B.3 first-order gradient methods that
can be used to compute solutions to Problem (27). The sequence θ m , defined via (27), leads
1
)-suboptimal solution in m many
to the solution of Problem (17) as m → ∞, providing a O( m
iterations if one uses standard proximal gradient descent methods; the convergence rate can be
improved to O( m12 ) if one uses the accelerated gradient descent version of the algorithm.
Instead of choosing f1 (θ) ≡ 0 one may also choose f1 (θ) = τ2 kθk22 , for a small value of τ > 0.
Interestingly, for small values of τ the minimizer to Problem (26) is also a minimizer of the
problem with the choice f1 (θ) = 0. This equivalence of solutions which holds true in much
more generality is often known as exact regularization of convex programs in the mathematical
programming literature — see for example [22]. Even if the two problems are not equivalent, the
choice of f1 (θ) = τ2 kθk22 always serves as an approximate solution to Problem (17). With this

39

choice of f1 (θ), one needs to solve a problem of the form:
p

min
θ

s.t.

X
τ
kθk22 +
wi |θi |
2
i=1

>

kX (y − Xθ)k∞ ≤ δ.

A solution to the above problem can be computed by considering its dual and applying (accelerated) proximal gradient methods on the dual formulation, as described in Section B.3. In this
approach, a two-stage iterative algorithm of the form (27) described above, is not required.
Algorithm 2 suggests that we solve Problem (17) repeatedly for different values of γ — it turns
out that the overall cost for solving all these problems is quite small. This is because (a) the
problems do not change much across different values of γ; and (b) for a fixed γ, while moving
across different values of k, the linear optimization problems are quite similar since the weights
|ρ0γ (|βik |)| do not change much across k. Thus the solutions obtained from one linear optimization
problem can be used as a warm-start to solve the next linear optimization problem. This is
found to reduce the overall computation time. Both the first-order methods (described above)
and simplex methods can gracefully take advantage of warm-starts.

B.3

Dual Gradient Method

Here we describe how to solve a problem of the form:
min
α

1
2 kα

−

ck22

+

p
X

wi |αi |

i=1

(28)

s.t. kAα − bk∞ ≤ δ,
where we assume that wi ≥ 0 and the set {α : kAα − bk∞ ≤ δ} is nonempty. Note that the
constraint set in (28) makes solving the primal form (28) challenging. However, due to the strong
convexity of the objective a dual is smooth (has Lipschitz continuous gradient) and the non
differentiability nicely separates across the dual variables – we thus use dual proximal gradient
algorithms to optimize (28). This trick is often used in optimization and signal processing, see
for example [17].
To derive a dual for Problem (28), we note that it can be written equivalently as:
min
α,ζ

1
2 kα

−

ck22

+

p
X
i=1

s.t. kζk∞ ≤ δ,
ζ = Aα − b.

40

wi |αi |

The minimum of the above problem can be obtained by maximizing a dual problem, obtained
by dualizing the equality constraints ζ = Aα − b; this consequently leads to the following
problem:

1
2
g(µ) :=
min
2 kα − ck2 + hµ, ζ − (Aα − b)i
ζ,α:kζk∞ ≤δ

p
P
+
wi |αi | .
i=1

The above can be simplified to:
g(µ)
= min
α

1
2 kα

−

ck22

− kµk1 δ − hµ, (Aα − b)i +

p
X

!
wi |αi |

i=1

=g1 (µ) − δkµk1 ,
where,
g1 (µ) = min
α

1
2 kα

− ck22 − hµ, (Aα − b)i +

p
X

!
wi |αi | .

i=1

b the unique minimizer of the above problem, is given by:
Note that α,
b = arg min
α
α

1
2 kα

− (c + A> µ)k22 +

p
X

wi |αi |,

i=1

i.e., α
bi = sgn(ci +

a>
i µ)

n
o
· max |ci + a>
µ|
−
w
,
0
,
i
i

for i = 1, . . . , p, where ai is the ith column of A. It follows from standard convex analysis [39]
that the function µ 7→ g1 (µ) is differentiable with its gradient given by:
b − b),
∇g1 (µ) = −(Aα
and its gradient is Lipschitz continuous:
k∇g1 (µ) − ∇g1 (µ0 )k ≤ kAk22 kµ − µ0 k,
where kAk2 denotes the largest singular value of A and for a vector u, the term kuk denotes the
usual `2 -norm of u.
By using standard quadratic programming duality theory [11], the minimum of Problem (28) can
be obtained by maximizing the unconstrained dual problem g(µ) in the dual variable µ, which
is equivalent to the following minimization problem:
min −g(µ) = min (−g1 (µ) + δkµk1 ) .
µ

µ

41

(29)

This problem is of the composite form [36], and proximal gradient descent methods [37, 36, 38]
apply to it directly.
For the special case of Problem (25), the method described above applies with wi = 0, i = 1, . . . , p.
Clearly, (29) is an `1 -regularized quadratic program, with the primal dual relationship being:
α = c + A> µ.
Note that Problem (28) needs to be solved several times during the course of Algorithm 1 and
Algorithm 2, across the different iterations. Fortunately, these problems are not completely
unrelated, in fact, they are quite “similar”. In Algorithm 2 the weights wi change; and in
Algorithm 1, the parameter c changes. Since the problems are similar, it is not unreasonable to
expect that the optimal dual variables corresponding to these two problems do not change much.
Thus, it is useful to initialize the dual variable µ for one instantiation of Problem (28) with the
(dual) solution obtained from another instantiation of Problem (28). This simple strategy leads
to substantial performance gains over solving the problems independent of one another.

B.4

Proof of Theorem 5

Proof. Note that the sequence β k , defined via (17) satisfies the following relationship:
h(β k ) =h(β k ; β k )
≥ min h(β; β k ) s.t. kX> (y − Xβ)k∞ ≤ δ
β

=h(β k+1 ; β k ).
Observing that
h(β k+1 ; β k ) ≥ h(β k+1 ),
we have:
h(β k ) = h(β k ; β k ) ≥ h(β k+1 ; β k ) ≥ h(β k+1 ),

(30)

and, thus, the sequence h(β k ) is decreasing. Subtracting h(β k ) from all sides of the above
inequality, we derive:
0 ≥ h(β k+1 ; β k ) − h(β k ) ≥ h(β k+1 ) − h(β k ).
The first part of the above display gives us, using (16):
0 ≥h(β k+1 ; β k ) − h(β k )
p D
E
X
=
ρ0γ (|βik |), |βik+1 | − |βik | = ∆(β k ),
i=1

which means ∆(β k ) ≤ 0 for all k. If ∆(β k ) < 0, then β k+1 leads to a strictly improved value

42

of the objective function. If ∆(β k ) = 0, then β k is a fixed point of the above update equation.
Hence, ∆(β k ) is a measure of how far β k is from a first-order stationary point of Problem (15).
The display in (30) shows that the objective values are decreasing, and, because the objective
values are all bounded below (by zero), the decreasing sequence converges.
In addition, we have that
h(β k ) − h(β k+1 ) ≥ −∆(β k ).
Adding the above for k = 1, . . . , K, we have:
h(β 1 ) − h(β K+1 ) ≥

P 
−∆(β k )
K≥k≥1

≥ K min

1≤k≤K

(31)


k

−∆(β ) ,

which leads to the following convergence rate:
n
o 1

−∆(β k ) ≤
h(β 1 ) − h(β K+1 )
1≤k≤K
K

1 
≤
h(β 1 ) − b
h ,
K
min

(32)
(33)

where (33) follows from (32) by using the observation that h(β k ) ↓ b
h.

B.5

Additional details on Algorithm 3

We seek an upper bound to a simple variant of Problem (15):
min h(β) :=
β

p
X

ργ (|βi |)

i=1
>

s.t. kX (y − Xβ)k∞ ≤ δ

(34)

βi = 0, i ∈ I c ,
b (1) ) := {i : βb(1) 6= 0, i = 1, . . . , p} ⊂ I, and I c is the complement of I. We assume,
where Supp(β
i
of course, that the feasible set in Problem (34) is nonempty. A simple method for constructing
I, which we found to be quite useful in practice, is presented below. Let B ⊂ {1, . . . , p} and B c
denote its complement. We define the following set:
n
o
F(B) := β : kX> (y − Xβ)k∞ ≤ δ, βi = 0, i ∈ B c .

43

b (1) be the solutions produced by Algorithm 1. Suppose we let B denote the support of
b (1) , β
Let α
b (1) ; the size of B is typically much smaller than p. If F(B) is nonempty, we take I = B. Note,
β
b (1) are only approximately equal: α
b (1) .
b (1) , β
b (1) ≈ β
however, that F(B) may be empty, because α
In this case, we need expand the set B, so that the set F(B) becomes nonempty. There may
be several ways to do this, but we found the following simple method to be quite useful in our
numerical experiments.
(1)

1. If F(B) is empty, we consider the set {|b
αi |, i ∈ B c } and find the index of the largest
(1)
element in this set, which we denote by: bi ∈ arg maxi∈B c |b
αi |.
2. Make B larger by including this new feature bi: we thus have B ← B ∪ {bi}.
3. Check if the resulting set F(B) is nonempty, if not, we repeat the above steps until F(B)
becomes nonempty.
4. We let I be the resulting set B obtained upon termination: I = B.
Problem (34), which is an optimization problem with fewer variables than Problem (15) is found
to deliver solutions that are better upper bounds to Problem (2). This also leads to better and
numerically more robust solutions than those available directly from Algorithm 1. The general algorithmic framework via sequential linear optimization, presented in Section 4.1.2, readily applies
to obtain good upper bounds to Problem (34).

B.6

Tighter bounds on β̂i ’s

The bounds described via (22) can be sharpened by making use of good upper bounds to the
solution of Problem (2). Towards this end, we need to reformulate (6). Note that in Problem (6),
if we take MU to be to be sufficiently large then this will lead to a solution for Problem (2). We
rewrite Problem (6) as follows:
min

α

β,z,α

s.t.

p
P

zi ≤ α

i=1

−δ ≤ dj − hqj , βi ≤ δ,

j = 1, . . . , p

−MU zj ≤ βj ≤ MU zj ,

j = 1, . . . , p

zj ∈ {0, 1},

j = 1, . . . , p,

where the optimization variables are β, z ∈ Rp and α ∈ R.

44

(35)

For a fixed α, consider the feasible set of Problem (35):




Sα = (β, z) :




p
X

>

zj ≤ α, kX (y − Xβ)k∞




≤ δ

j=1

|βj | ≤ MU zj , zj ∈ {0, 1}, j = 1, . . . , p

.





Observe that
Sα ⊂ S α ,
where





S α = (β, z) :




p
X

(36)

zj ≤ α, kX> (y − Xβ)k∞




≤ δ

j=1

|βj | ≤ MU zj , zj ∈ [0, 1], j = 1, . . . , p





is obtained by relaxing the binary variables zj ∈ {0, 1} into the continuous variables zj ∈ [0, 1], for
all j = 1, . . . , p. Noting that Sα ⊂ Sα0 for α ≤ α0 ; and using this along with (36) we have
Sα∗ ⊂ Sα0 ⊂ S α0 ,
where α∗ is the optimum value, and α0 is an upper bound to Problem (35), and hence α0 ≥ α∗ .
The above inequality leads to the following chain of inequalities:
min
(β,z)∈Sα∗

max
(β,z)∈Sα∗

βi ≥
βi ≤

min
(β,z)∈Sα0

max
(β,z)∈Sα0

βi ≥

min
(β,z)∈S α0

βi ≤

max
(β,z)∈S α0

βi := µ−
i (α0 )

βi := µ+
i (α0 ).

+
The quantities at the right end above, i.e. µ−
i (α0 ) and µi (α0 ), can be computed by solving a
pair of linear optimization problems:

µ+
i (α0 ) :=
s.t.

max βi
β

kX> (y − Xβ)k∞ ≤ δ,
kβk∞ ≤ MU ,
kβk1 ≤ MU α0 ,
(37)

µ−
i (α0 ) :=
s.t.

min βi
β

kX> (y − Xβ)k∞ ≤ δ,
kβk∞ ≤ MU ,
kβk1 ≤ MU α0 .

45

+
b
The quantities µ−
i (α0 ) and µi (α0 ) are lower and upper bounds, respectively, for βi — the bounds
 +
−
depend upon α0 and MU . Note that µi (α0 ) := max µi (α0 ), −µi (α0 ) provides an upper bound
b ∞ — this suggests a way to
to |βbi |, which consequently leads to an improved estimate for kβk
−
+
adaptively refine MU , and, thus, µi (α0 ) and µi (α0 ).

C

Additional Experiments

This section complements the experimental results shown in the main body of the paper. Table 4
is an elaborate version of the representative results displayed in Figure 5. Here, we consider
different values of SNR and also display the prediction errors. The results show that Discrete
Dantzig Selector outperforms the `1 -Dantzig Selector based methods in terms of estimating the
true underlying regression coefficients, and does so with better variable selection properties.
Example-A (n = 200,p = 500)
Metric
kβb − β∗k22
Variable Selection Error
Prediction Error
Number of Nonzeros
kβb − β∗k22
Variable Selection Error
Prediction Error
Number of Nonzeros

SNR
3
3
3
3
10
10
10
10

Metric
kβb − β∗k22
Variable Selection Error
Prediction Error
Number of Nonzeros
kβb − β∗k22
Variable Selection Error
Prediction Error
Number of Nonzeros

SNR
3
3
3
3
10
10
10
10

L0-DS
3.858 (0.466)
6 (0.543)
0.158 (0.013)
24 (0.274)
0.635 (0.053)
1 (0.137)
0.031 (0.002)
21 (0.137)

L0-DS-Pol
1.123 (0.254)
16 (0.844)
0.581 (0.038)
35 (0.663)
0.172 (0)
0 (0)
0.249 (0)
20 (0)

Example-B (n = 200,p = 500)
L1-DS
5.266 (0.162)
102 (1.327)
0.18 (0.004)
122 (1.327)
1.58 (0.049)
102 (1.327)
0.054 (0.001)
122 (1.327)

Warm
3.429 (0.19)
43 (1.327)
0.133 (0.004)
63 (1.266)
0.664 (0.035)
29 (1.371)
0.028 (0.001)
49 (1.371)

L1-DS
6.452 (0.263)
43 (1.538)
0.172 (0.004)
47 (1.658)
3.245 (0.097)
38 (1.357)
0.056 (0.002)
48 (1.357)

Warm
7.674 (0.347)
21 (1.116)
0.193 (0.008)
22 (1.096)
2.529 (0.357)
27 (1.343)
0.051 (0.004)
35 (1.266)

Metric
kβb − β∗k22
Variable Selection Error
Prediction Error
Number of Nonzeros
kβb − β∗k22
Variable Selection Error
Prediction Error
Number of Nonzeros

SNR
3
3
3
3
10
10
10
10

Example-C (n = 100,p = 500)
L0-DS
7.322 (0.65)
10 (0.573)
0.168 (0.015)
13 (0.271)
2.395 (0.379)
10 (0.795)
0.054 (0.004)
18 (0.685)

L0-DS-Pol
6.823 (0.693)
9 (0.814)
0.966 (0.051)
13 (0.392)
1.988 (0.333)
10 (0.795)
0.515 (0.036)
18 (0.685)

L0-DS
150.386 (8.993)
11 (0.247)
0.155 (0.007)
16 (0.362)
45.599 (2.508)
9 (0.411)
0.045 (0.001)
17 (0.151)

L0-DS-Pol
64.17 (4.888)
16 (0.573)
0.109 (0.003)
23 (0.362)
13.773 (1.39)
11 (0.482)
0.052 (0.001)
23 (0.392)

L1-DS
137.073 (4.511)
38 (1.266)
0.132 (0.005)
46 (1.447)
47.113 (1.016)
36 (1.387)
0.043 (0.001)
50 (1.266)

Warm
102.583 (5.711)
14 (0.658)
0.123 (0.005)
22 (0.693)
38.256 (2.068)
14 (0.724)
0.035 (0.001)
26 (0.663)

Example-D (n = 100,p = 300)
Metric
kβb − β∗k22
Variable Selection Error
Prediction Error
Number of Nonzeros
kβb − β∗k22
Variable Selection Error
Prediction Error
Number of Nonzeros

SNR
3
3
3
3
10
10
10
10

L0-DS
L0-DS-Pol
L1-DS
Warm
0.019 (0.006) 0.015 (0) 0.191 (0.027) 0.018 (0.007)
0 (0)
0 (0)
36 (0.693) 12 (0.85)
0.033 (0.005) 0.955 (0) 0.101 (0.007) 0.028 (0.003)
2 (0)
2 (0)
38 (0.693) 14 (0.85)
0.005 (0.002) 0.004 (0) 0.057 (0.008) 0.006 (0.002)
0 (0)
1 (0.082) 36 (0.693) 13 (0.877)
0.006 (0.002) 0.527 (0.006) 0.03 (0.002) 0.007 (0.001)
2 (0)
3 (0.082) 38 (0.693) 15 (0.877)

Table 4: Tables showing the statistical performance of four different methods: “L0-DS”, “L0-DS-Pol”,
“L1-DS” and “Warm”, described in Section 7. The standard errors are given in parentheses (they are
computed in the same fashion as in Figure 5). The Discrete Dantzig Selector based methods deliver
models with good accuracy in estimating the regression coefficients, and the estimated models are sparser
than those for the `1 -based method and the method based on Algorithm 2, which is a heuristic strategy
to approximate good upper bounds for the Discrete Dantzig Selector problem.

An important advantage of the Discrete Dantzig Selector based methods is that they deliver
models that are very sparse. The polished version of the Discrete Dantzig Selector is found
to exhibit better statistical performance than the original Discrete Dantzig Selector estimator.
Table 5 compares the polished versions of the Discrete Dantzig Selector and the `1 -Dantzig
Selector, and finds that the performance of the former approach is significantly better.
46

Example-A (n = 200, p = 500)
Metric
b
kβ − β ∗ k22
Variable Selection Error
Prediction Error
Number of Nonzeros
b − β ∗ k2
kβ
2
Variable Selection Error
Prediction Error
Number of Nonzeros

SNR
3
3
3
3
10
10
10
10

L0-DS-Pol
1.123 (0.254)
16 (0.844)
0.581 (0.038)
35 (0.663)
0.172 (0)
0 (0)
0.249 (0)
20 (0)

L1-DS-Pol
2.163 (0.139)
30 (1.116)
0.777 (0.02)
50 (1.146)
0.282 (0.011)
11 (0.392)
0.292 (0.003)
31 (0.392)

Table 5: Tables comparing the polished version of Discrete Dantzig Selector with the polished version
of `1 -Dantzig Selector. The standard errors are given in parentheses. The statistical performance of
`1 -Dantzig Selector is inferior, most likely due to its weaker variable selection properties.

47

References
[1] Top500 Supercomputer Sites, Directory page for Top500 lists. Result for each list since June 1993.
http://www.top500.org/statistics/sublist/. Accessed: 2013-12-04.
[2] A. Beck and M. Teboulle. A fast iterative shrinkage-thresholding algorithm for linear inverse problems.
SIAM Journal on Imaging Sciences, 2(1):183–202, 2009.
[3] S. R. Becker, E. J. Candès, and M. C. Grant. Templates for convex cone problems with applications
to sparse signal recovery. Mathematical Programming Computation, 3(3):165–218, 2011.
[4] D. P. Bertsekas. Nonlinear Programming. Athena Scientific, Belmont, Massachusetts, 2nd edition,
1999.
[5] D. Bertsimas and A. King. Or forum – an algorithmic approach to linear regression. Operations
Research, 2015.
[6] D. Bertsimas and R. Mazumder. Least quantile regression via modern optimization. Annals of
Statistics, 42(6):2494–2525, 2014.
[7] D. Bertsimas and R. Weismantel. Optimization over integers. Dynamic Ideas Belmont, 2005.
[8] D. Bertsimas, A. King, and R. Mazumder. Best subset selection via a modern optimization lens.
Annals of Statistics, 44(2):813–852, 2016.
[9] P. Bickel, Y. Ritov, and A. Tsybakov. Simultaneous analysis of lasso and dantzig selector. Annals of
Statistics, 37:1705–1732, 2009.
[10] R. E. Bixby. A brief history of linear and mixed-integer programming computation. Documenta
Mathematica, Extra Volume: Optimization Stories, pages 107–121, 2012.
[11] S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University Press, Cambridge, 2004.
[12] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein. Foundations and Trends in Machine Learning.
Number 3(1). Now Publishers, 2011.
[13] P. Bühlmann and S. van-de-Geer. Statistics for high-dimensional data. Springer, 2011.
[14] S. Burer and A. Saxena. The MILP road to MIQCP. In Mixed Integer Nonlinear Programming, pages
373–405. Springer, 2012.
[15] E. Candes, M. Wakin, and S. Boyd. Enhancing sparsity by reweighted `1 minimization. Journal of
Fourier Analysis and Applications, 14(5):877–905, 2008.
[16] E. Candes and T. Tao. The Dantzig selector: statistical estimation when p is much larger than n.
Annals of Statistics, pages 2313–2351, 2007.
[17] P. L. Combettes, DJ. Dũng, and B. C. Vũ. Dualization of signal recovery problems. Set-Valued and
Variational Analysis, 18(3-4):373–404, 2010.
[18] D. Donoho and I. Johnstone. Ideal spatial adaptation by wavelet shrinkage. Biometrika, 81:425–455,
1994.
[19] B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani. Least angle regression (with discussion). Annals
of Statistics, 32(2):407–499, 2004. ISSN 0090-5364.
[20] J. Fan and R. Li. Variable selection via nonconcave penalized likelihood and its oracle properties.
Journal of the American Statistical Association, 96(456):1348–1360(13), 2001.
[21] R. M. Freund, P. Grigas, and R. Mazumder. A new perspective on boosting in linear regression via
subgradient optimization and relatives. Annals of Statistics (to appear), 2017.
[22] M. P. Friedlander and P. Tseng. Exact regularization of convex programs. SIAM Journal on Optimization, 18(4):1326–1350, 2007.

48

[23] J. Friedman, T. Hastie, H. Hoefling, and R. Tibshirani. Pathwise coordinate optimization. Annals of
Applied Statistics, 2(1):302–332, 2007.
[24] F. Giannessi and E. Tomasin. Nonconvex quadratic programs, linear complementarity problems,
and integer linear programs. In 5th Conference on Optimization Techniques Part I, pages 437–449.
Springer, 1973.
[25] I. Gurobi Optimization. Gurobi optimizer reference manual, 2015. URL http://www.gurobi.com.
[26] T. Hastie, R. Tibshirani, and J. Friedman. The Elements of Statistical Learning, Second Edition:
Data Mining, Inference, and Prediction (Springer Series in Statistics). Springer New York, 2 edition,
2009. ISBN 0387848576.
[27] T. Hastie, R. Tibshirani, and M. Wainwright. Statistical Learning with Sparsity: The Lasso and
Generalizations. CRC Press, FL, 2015.
[28] R. Hemmecke, M. Köppe, J. Lee, and R. Weismantel. Nonlinear integer programming. In 50 Years
of Integer Programming 1958-2008, pages 561–618. Springer, 2010.
[29] G. M. James and P. Radchenko. A generalized Dantzig selector with shrinkage tuning. Biometrika,
96:323–337, 2009.
[30] G. M. James, P. Radchenko, and J. Lv. Dasso: connections between the Dantzig selector and lasso.
Journal of the Royal Statistical Society: Series B (Statistical Methodology), 71(1):127–142, 2009.
[31] M. Jünger, T. M. Liebling, D. Naddef, G. L. Nemhauser, W. R. Pulleyblank, G. Reinelt, G. Rinaldi,
and L. A. Wolsey. 50 Years of Integer Programming 1958-2008: From the Early Years to the Stateof-the-art. Springer Science & Business Media, 2009.
[32] J. T. Linderoth and A. Lodi. MILP software. Wiley encyclopedia of operations research and management science, 2010.
[33] H. Liu, T. Yao, and R. Li. Global solutions to folded concave penalized nonconvex learning. Annals
of Statistics, 44 (2):629–659, 2016.
[34] Z. Lu and Y. Zhang. Sparse approximation via penalty decomposition methods. SIAM Journal on
Optimization, 23(4):2448–2478, 2013.
[35] R. Mazumder, J. Friedman, and T. Hastie. Sparsenet: Coordinate descent with non-convex penalties.
Journal of the American Statistical Association, 117(495):1125–1138, 2011.
[36] Y. Nesterov. Gradient methods for minimizing composite functions. Mathematical Programming, 140
(1):125–161, 2013.
[37] Y. Nesterov. Introductory Lectures on Convex Optimization: A Basic Course. Kluwer, Norwell, 2004.
[38] N. Parikh and S. Boyd. Proximal algorithms. Foundations and Trends in optimization, 1(3):123–231,
2013.
[39] R. Rockafellar. Convex Analysis. Princeton University Press, Princeton, 1996.
[40] R. Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical
Society, Series B, 58:267–288, 1996.
[41] J. P. Vielma. Mixed integer linear programming formulation techniques. SIAM Review, 57(1):3–57,
2015.
[42] H. P. Williams. Model building in mathematical programming. John Wiley & Sons, 2013.
[43] C.-H. Zhang and J. Huang. The sparsity and bias of the lasso selection in high-dimensional linear
regression. Annals of Statistics, 36(4):1567–1594, 2008.

49

