arXiv:1609.05609v1 [math.ST] 19 Sep 2016

Selective sampling after solving a
convex problem
Xiaoying Tian Harris, Snigdha Panigrahi, Jelena Markovic, Nan Bi,
Jonathan Taylor
Abstract: We consider the problem of selective inference after solving a
(randomized) convex statistical learning program in the form of a penalized
or constrained loss function. Our first main result is a change-of-measure
formula that describes many conditional sampling problems of interest in
selective inference. Our approach is model-agnostic in the sense that users
may provide their own statistical model for inference, we simply provide
the modification of each distribution in the model after the selection.
Our second main result describes the geometric structure in the Jacobian
appearing in the change of measure, drawing connections to curvature measures appearing in Weyl-Steiner volume-of-tubes formulae. This Jacobian
is necessary for problems in which the convex penalty is not polyhedral,
with the prototypical example being group LASSO or the nuclear norm.
We derive explicit formulae for the Jacobian of the group LASSO.
To illustrate the generality of our method, we consider many examples
throughout, varying both the penalty or constraint in the statistical learning problem as well as the loss function, also considering selective inference
after solving multiple statistical learning programs. Penalties considered include LASSO, forward stepwise, stagewise algorithms, marginal screening
and generalized LASSO. Loss functions considered include squared-error,
logistic, and log-det for covariance matrix estimation.
Having described the appropriate distribution we wish to sample from
through our first two results, we outline a framework for sampling using
a projected Langevin sampler in the (commonly occuring) case that the
distribution is log-concave.

1. Introduction
Based on the explosion of freely available and high quality statistical software,
practicing data scientists or statisticians can rely on untold numbers of statistical learning methods that explore their data. Having found an interesting
pattern through such methods, in order to report their findings in the scientific literature, the data scientist is confronted in defending the significance of
their findings. Naively assessing significance using methods that ignore their
exploration is recognized as flawed. A common, though not often used, solution
is to use data splitting to evaluate the significance. There has been significant
recent research on developing methods of inference that have the same type of
guarantees as data splitting but offer more power. These methods are generally
referred to as methods for selective inference (Benjamini, 2010; Berk et al., 2013;
Barber and Candes, 2014; Fithian et al., 2014; Lee et al., 2016; Tian et al., 2015;
Tian and Taylor, 2015). Loosely speaking, selective inference recognizes the in1
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

2

herent selection biases in reporting the most “significant” results from various
statistical models and attempts to adjust for the bias in a rigorous framework.
At a high level, selective inference involves two stages: first, query the data
by applying some function. Often this function might be the solution of some
convex optimization problem. Second, posit a model based on the outcome of
the query and perform inference for the parameters or statistical functionals in
such models. In order to perform valid inference in this second stage, a common approach espoused above is to condition on the result of the query. An
alternative approach is to address the problem via reduction to a problem of simultaneous inference Berk et al. (2013). Under the conditional approach, when
the query involves solving a convex problem, we are interested in distributions
of the form
S|β̂(S) ∈ A
(1)
where S ∼ F represents our data, β̂ denotes the solution to a convex optimization problem and A is some event of interest, set by the data analyst.
In this paper, we address computational problems in the conditional approach
to selective inference. In adopting this conditional approach, it quickly becomes
apparent distributions of the form (1) can be rather complex. In this work,
we propose a unified sampling approach that is feasible for a wide variety of
problems. Combined with the randomization idea in Tian and Taylor (2015),
this significantly improve both the applicability of selective inference in practice,
and the power of the selective tests.
1.1. A canonical example
As a concrete example, we consider solving the LASSO Tibshirani (1996a)
at some fixed λ, as considered in Lee et al. (2016). Suppose we observe data
(X, y), with X ∈ Rn×p and y ∈ Rn . The query by the data analyst returns
sign(β̂(X, y; λ)) (with sign(0) = 0) where
1
β̂(X, y; λ) = argmin ky − Xβk22 + λkβk1 .
2
β

(2)

The sign vector can be expressed in terms of Ê, the support of β̂ and ẑE the
signs of the non-zero coefficients of β̂. In principle, after observing (Ê, ẑE ) =
(Eobs , zE,obs ), there is no restriction on the choice of models. For simplicity,
we posit a parametric normal model that is based on E (dropping the ‘obs′
notation) and fixed X,

2
ME = N (XE bE , σE
I) : bE ∈ RE
(3)

2
with σE
known. Target of inference is naturally bj|E , the parameters of this
model. Selective inference essentially takes any distribution F ∈ ME and consider its corresponding selective distribution by conditioning on {(Ê, ẑE ) =
(Eobs , zE,obs )}. In this case, our selective distribution is exactly of the form

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

3

of (1), with A being the quandrant specified by the nonzero coefficients E and
their signs zE .
Despite the simple form of A, the constraint induced a much more intricate
set in the space of y. In this special case, Lee et al. (2016) has worked out the
constraint set on y and selective tests can be computed explicitly. However, this
is not generally true, especially in the context of randomized selection procedures
proposed by Tian and Taylor (2015); Fithian et al. (2014). In many cases, the
formsof selective distributions are not explicitly computable and we have to use
sampling methods to approximate the conditional distribution. Given a large
enough sample from such distribution, we can construct the selective tests. In the
following paragraph, we first introduce the randomized Lasso procedure which is
a special case of randomized selection procedures proposed in Tian and Taylor
(2015). Then we would illustrate what is the appropriate conditional distribution
to sample and given a large sample from this distribution, how we carry out the
tests for bj|E .
First, we propose to incorporate randomness into (2). In particular, we consider solving
ǫ
1
β̂(X, y, ω; λ) = argmin ky − Xβk22 + λkβk1 − ω T β + kβk22
2
2
β

(4)

where ω ∼ G is a random vector independent of (X, y), whose distribution is
chosen by the data analyst, hence known. We will assume that G is supported
on all of Rp with density g. The ridge term with small parameter ǫ ensures the
problem above has a solution. This resembles the ridge term in the elastic net
proposed by Zou and Hastie (2005).
This type of randomized convex program has been considered in Tian and Taylor
(2015). It is espoused as it significantly increase the power in the inference stage
without much loss in the quality of the selected model. In particular, increasing
the scale of the randomization inevitably deteriorates the model selection quality, similar to using a smaller training set to choose a model in the context of
data splitting. However, even a small amount of randomization has empirically
shown a fairly noticable increase in power. Another benefit of randomization
include a form of robustness to rare selection events in the selective CLT of
Tian and Taylor (2015).
After solving (4), we use its solution β̂(X, y, ω; λ) for model selection. Similar
to unrandomized version, the appropriate law after model selection is
y | β̂(X, y, ω; λ) ∈ A,

y ∼ F,

where F is a member of the parametric normal model specified by (3).
Various selective tests can be constructed based on this law including the
goodness-of-fit tests and the p-values and confidence intervals for bj|E s. We
discuss how to construct tests for bj|E s here as they are a special case of testing a
single parameter in a multi-parameter exponential family. The general approach
is laid out in Fithian et al. (2014).

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

4

Selective inference recognizes the fact that interest in bj|E is due to the query
returning E as the nonzero variables of the solution to (4). In order to provide
valid inference, we therefore consider only the part of sample space that will yield
the same model, i.e. conditioning on {(Ê, ẑE ) = (Eobs , zE,obs )}. Conditioning
each model in ME yields a new parametric model, the selective model


dF ∗
(·) ∝ 1{(Ê,ẑE )=(Eobs ,zE,obs ),F ∈ME } (·) .
(5)
M∗(Eobs ,zE,obs ) = F ∗ :
dF
where (Eobs , zE,obs ) are the observed variables and signs.
Note that this approach is not tied to the parametric model in (3). In particular, when the parametric normal model is not appropriate, one might replace
M with some other statistical model (i.e. some other collection of distributions)
and carry out statistical inference in this new selective model. This approach
was laid out in detail in Fithian et al. (2014), with some asymptotic justification
for nonparametric models developed in Tian and Taylor (2015). In fact, in the
Lee et al. (2016) the (pre-conditioning) statistical model proposed was not ME ,
rather it was the saturated model

(6)
M = N (µ, σ 2 I) : µ ∈ Rn .

A fairly simple calculation shows that the hypothesis tests and confidence intervals are exactly the same whether we had started with (6) or (3). However
this is not the case for examples like forward stepwise, discussed below and in
Fithian et al. (2015); Tibshirani et al. (2014).
Let us now consider how we might go about inference in (5). For instance,
suppose we are simply interested in a goodness-of-fit test and we will consider
the model M = {F0 } so that M∗ consists of F0 restricted to the event that fixes
the active set and signs of the LASSO to be (Eobs , zE,obs ). Call this distribution
F0∗ . A natural way to test H0 : F = F0 would be follow Fisher’s approach by
choosing some test statistic T = T (X, y) and compare our observed value Tobs
to the distribution of T under F0∗ . In order to do this, it is sufficient to describe
the law β̂λ,∗ (F0 ), i.e. the push forward of F0 and condition this law on the event
that fixes the active set and signs to be (Eobs , zE,obs ). Unfortunately, this is no
easy task. Of course, a natural alternative is to use Monte Carlo. A large enough
sample from F0∗ is sufficient to carry out this goodness-of-fit test. A version of
this goodness-of-fit test is considered in Section 4.7 below.
The issue of what Monte Carlo method to use remains. Our main contribution
in this work is an explicit description of how to sample from distributions such
as F0∗ . To be precise, we describe how to sample from conditional distributions
where the conditioning depends on the solution to a convex problem. A naive
way to sample from this distribution would be to take a Monte Carlo sampler to
draw from F0 , retaining only those points where the active set and signs of the
LASSO agree with (Eobs , zE,obs ). This requires solving (2) at each sample point,
and is clearly infeasible. We might say that this is a direct way to condition the
push-forward distribution β̂λ,∗ (F0 ) on the selection event.

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

5

Our approach here relies on what is essentially the pull-back of the conditional
distribution itself. In concrete terms, we construct an explicit inverse map to
problems like (2) on the selection event. Then, instead of sampling (1), we
realize the law F0∗ by sampling β̂λ together with some auxillary variables on the
conditioning event. After aquiring these samples, we can use the inverse map
to reconstruct y that is essentially distributed according to the law F0∗ . The
construction for (2) is described in Section 2.2.
We will see that the pull-back construction is particularly simple for the randomized lasso (4). The randomization was inspired by the differential privacy
literature Dwork et al. (2015), though we will see below that it is already very
similar to data splitting. Besides making the pull-back simple to compute, we advocate the use of randomization in the query stage in that conditional inference
in the second stage is often more powerful if the selection stage is carried out
with randomization than without Fithian et al. (2014); Tian and Taylor (2015).
Increasing the scale of the randomization inevitably deteriorates the model selection quality, similar to using a smaller training set to choose a model in the
context of data splitting. However, even a small amount of randomization has
empirically shown a fairly noticable increase in power. Another benefits of randomization include a form of robustness to rare selection events in the selective
CLT of (Tian and Taylor, 2015).
Inference for linear functionals is slightly more complicated than the simple
goodness-of-fit test. To conclude this section, let us describe the form of the
pull-back construction when we have solved (4) and observed the active set and
signs (E, zE ) (dropping the ‘obs′ notation). In principle, there is no restriction
on the choice of models. For simplicity, we posit a parametric normal model
that is based on a subset of variables Ē:
o
n
2
MĒ = N (XĒ βĒ , σĒ
I) : βĒ ∈ RĒ

2
with σĒ
known. Often, the data analyst may opt to take Ē = Eobs , though
this is not strictly necessary. Based on (E, zE ) she might consult the relevant
literature (or data independent of y) and choose to include or delete some variables from the set E. When the variance is unknown, one can either plug in
2
a consistent estimate of σE
(c.f. Lemma 14 of Tian and Taylor (2015)) or in2
clude the parameter σĒ in the model and solve instead the square-root Lasso
Belloni et al. (2014). If a parametric model is not appropriate then one might
use a normal approximation for the pair X T y, and rewrite the selection event
in terms of this statistic as opposed to just y. Such an approach is considered
in Tian and Taylor (2015), appealing to the selective CLT.
Having fixed model MĒ , and knowing the density g for the randomization
ω, suppose the data analyst now wants to test H0 : βj|Ē under the assumption
F ∈ MĒ . Let PĒ\j denote orthogonal projection onto col(XĒ\j ). In the normal
z-test for βj|Ē , we essentially consider the distribution of

rj = y − µĒ\j ,

µĒ\j = PĒ\j y.

rj can be recognized as the residual from the model with variables Ē \ j. Selecimsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

6

tive tests for βj|Ē is essentially considering the law of rj under the conditional
distributions. An application of our main result Theorem 1 implies that this law
(with the appropriate augmentation variables) has density proportional to

 

 

1
ẑE
β̂E
T
2
+ǫ·
(rj , β̂E , ẑ−E ) 7→ exp − 2 kµĒ\j + rj k2 ·g X XE β̂E + λ ·
ẑ−E
0
2σĒ
(7)
supported on
n
o
col(PĒ\j )⊥ × (β̂E , ẑ−E ) : sign(β̂E ) = ẑE , kẑ−E k∞ ≤ 1 .
The augmentation variables (β̂E , ẑ−E ) are variables involved in the optimization
problem (4). Their support ensures we stay on the selection event of interest.
In fact, their support is recognizable as a subset of
n
o
(β̂, ẑ) : ẑ ∈ ∂(k · k1 )(β̂) ,
with ∂ denoting the subdifferential with respect to β̂. In particular, this is the
subset for which the non-zero coefficients of β̂ are E with signs ẑE .
In this simple case, the construction of the inverse map is remarkably simple.
If we write out the KKT condition for the optimization problem (4), we have
X T X β̂ − X T y + λ · ẑ − ω + ǫ · β̂ = 0
which can be rewritten as
ω = X T X β̂ − X T y + λẑ + ǫβ̂.
where (β̂, ẑ) are the optimization variables and the subgradient of (4). In this
sense, we have inverted the KKT conditions on the selection event of interest,
explicitly parameterizing the pairs (y, ω) that yield active set and signs (E, zE ).
We call this density the pull-back as this is recognizable precisely as the differential geometric pull-back of a differential form on the interior of the selection
event, a subset of the manifold Rn × Rp . As g is assumed to have a density, the
boundary of the selection event is ignorable. To those readers who do not recognize the term pull-back, Theorem 1 essentially just applies a standard change of
variables to realize a conditional law for (y, ω) via a law for (rj , βE , ẑ−E ). Note
that the unselective law of rj is that of y projected onto the residual space form
model Ē \ j.
The reader may ask what we have gained through this exercise. By constructing an inverse map for (y, ω) on the selection event, we do not have to check
whether the pair (y, ω) = (y, X T X β̂ − X T y + λ · ẑ + ǫ · β̂) satisfy the KKT
conditions for the pair (E, zE ). That is, we have removed the difficulty of conditioning from the problem. We are still left with a sampling problem, though we
have an explicit density as well as a fairly simple support. Further, the density
(7) is log-concave in (rj , βE , z−E ) with a simple support hence methods such as
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

7

projected Langevin Bubeck et al. (2015) may be run efficiently as the projection
operator is cheap to compute. Section 7 discusses several other examples.
Finally, and most importantly, the reader who may also be the data scientist
may ask what to do with samples from this density in order to carry out a
test of H0 : βj|Ē (F ) = 0. Let β̂j|Ē : Rn → R denote the map that computes
OLS coefficient j in the model with variables Ē. Then, given a sufficiently large
sample (rj,b , βE,b , z−E,b )B
b=1 from density (7), the data analyst will compare the
empirical distribution of

B
β̂j|Ē (µĒ\j + rj,b )

b=1


B
= β̂j|Ē (rj,b )

b=1

to the observed value β̂j|Ē (y). Selective confidence intervals can be constructed
by tilting the empirical distribution, though if the true parameter is far from
0, then a reference distribution other than (7) is perhaps more appropriate.
One might try replacing rj,⊥ in the density above with rj,⊥ − β̄j|E where β̄j|E
is an approximate selective MLE or pseudo MLE Panigrahi et al. (2016). We
do not pursue this further here, leaving this for future work. A data analyst
concerned about slow mixing in a given MCMC scheme to draw from (7) may
take some reversible MCMC algorithm and carry out the exact tests described
in Besag and Clifford (1989).
1.2. General approach
Having described what might be the canonical example, the LASSO with a
parametric Gaussian model for inference, we now lay out our general approach.
We consider a randomized version of the optimization problem in (1) formulated as
ǫ
β̂(S, ω) = argmin ℓ(β; S) + P(β) − ω T β + kβk22
(8)
p
2
β∈R
where ℓ is some smooth loss involving the data, P is some structure inducing
convex function, ǫ > 0 is some small parameter that is sometimes necessary in
order to assure the program has a solution and ω ∼ G is a randomization chosen
by the data analyst. Our main goal is to sample
(S, ω)|β̂(S, ω) ∈ A

(9)

Our reasons for considering the problem (8) rather than an unrandomized problem are described in the LASSO example above and also hold for the optimization and sampling problems.
Namely, we expect an increase in power in the second stage following even
a small randomization. Further, as in the LASSO example above, this randomization often allows us to cast the sampling problem as sampling from a
distribution on a space that is much simpler than if we had not randomized.
That is, the law induced by the pull-back measure in a randomized program is
often supported on a simpler region, as opposed to the push-forward measure
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

8

of the non-randomized program. For inference in the non-randomized case, one
might take the approach of sending the scale of randomization to 0, though we
do not pursue this here.
Another feature of the sampling problem related to (8) is that the sampler
somewhat decouples the statistical model from the optimization variables. In
this sense, our main result provides ways to sample in a model-agnostic fashion:
data analysts can supply their own model M, resulting in selective model M∗ .
Of course, for inference in M∗ , the analyst may have to use other techniques to
reduce their problem to sampling from a particular distribution in M∗ . In the
example above, standard exponential family techniques were used to eliminate
nuisance parameters βĒ\j|Ē (F ) and the problem was reduced to sampling from
only one distribution, constructed by conditioning a distribution in M∗ on the
sufficient statistic corresponding to the nuisance parameters.
By decoupling, we do not mean statistical independence in any sense. We
mean that each distribution in the corresponding selective model is supported
on Ω × C where Ω is the original probability space for our data S and C is a set
of optimization variables related to the structure inducing function P. Formally,
we should note that in the pull-back variables, the selective model describes the
distribution of tuples (s, β, z) rather than the original probability space Ω × Rp .
Hence, our sampler produces tuples (s, β, z) rather than pairs (s, ω) though ω
can always be reconstructed via the map ω = ∇ℓ(β; s) + z + ǫ · β.
The two most common examples of interest in statistical learning are
P(β) = hK (β) = sup ν T β
ν∈K

for some convex K ∋ 0, i.e. a seminorm. In this case, C is typically a subset of
the normal bundle of K (c.f. Adler and Taylor (2007); Schneider (1993) )
{(β, z) : β ∈ Nz K}
where Nz K is the normal cone of K at z, polar to the support cone of K at z.
The other common example is a constraint on a seminorm, i.e.
(
0 β∈K
P(β) = IK (β) =
(10)
∞ β 6∈ K
where K = {b : kbk ≤ 1} for some seminorm k·k. In this case C is again typically
a subset of N (K) of the form
{(β, z) : z ∈ Nβ K} .

(11)

Such sets arise naturally from the KKT conditions of (8):
ω = ∇ℓ(β; S) + z + ǫ · β.
Queries related to active sets of variables in the case of the LASSO or group
LASSO or the rank of matrices in the case of the nuclear norm correspond
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

9

to smooth subsets of the corresponding normal bundles. When the structure
inducing penalty is polyhedral in nature, the sets C are typically polyhedral.
For norms with curved unit balls this is no longer the case, and curvature comes
in to play. We treat the group LASSO as a canonical example of this in Section
5.
1.3. Related work
Most of the theoretical work on high-dimensional data focuses on consistency,
either the consistency of solutions Negahban et al. (2010); Van de Geer (2008)
or the consistency of the models Wainwright (2009); Zhao and Yu (2006).
In the post selection literature, Berk et al. (2013) proposed the PoSI approach, which reduce the problem to a simultaneous inference problem. Because
of the simultaneity, it prevents data snooping from any selection procedure, but
also results in more conservative inference. In addition, the PoSI method has
extremely high computational cost, and is only applicable when the dimension
p < 30 or for very sparse models. The authors Meinshausen et al. (2012) proposed a method for computing p-values that controls false discovery rate (FDR)
among all variables. The knockoff filter of Barber and Candes (2014) provides
similar control of FDR for all variables in the full model. What distinguishes the
conditional approach from these simulataneous approaches is that the hypotheses tested, or parameters for which intervals are formed, in selective inference
are chosen as a function of the data. Hence, the methods of inference are not
always directly comparable.
1.4. Outline of paper
We propose the main Theorem 1 of this paper in Section 3, putting forth the
sampling density conditional on a selection event. The highlight of this theorem is the reparametrization map that allows to reconstruct the randomization
as a function of optimization variables and data and allows us to compute the
required conditional density explicitly, up to a normalizing constant. This is followed by a variety of examples of convex optimization programs. The support
of the selective density in the first set of examples- LASSO and variants in Section B, graphical models in Section 4.8, forward stepwise in Section 4.4 can be
described by polyhedral geometry. We follow this up with a section describing
more complex problems where the Jacobian involves a curvature component, the
group lasso illustrated as the prototypical example. We extend the selective sampler to the case in which the data analyst considers multiple views, or queries,
of the data in Section 6. A selective version of Fisher’s exact test in which the
data analyst chooses the sufficient statistics based on the data is described in
Section 4.7. Finally, we advocate the projected Langevin sampling technique
Bubeck et al. (2015) in Section 7 to sample from a log-concave selective density as described in earlier examples. Each update in such an implementation
involves a projection onto a set of constraints induced by the selection event.
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

10

The computational cost of each step is often minimal due to a much simpler
constraint region using our reparametrization.
2. Inverting the optimization map
In this section, we consider the general problem of constructing an explicit
inverse to the solution of a convex problem. In this section, we focus on convex
problems without additional randomization, deferring randomization to Section
3.
Consider a statistical learning problem of the form
minimize
ℓ(β; S) + P(β),
p
β∈R

S ∼ F,

(12)

where F is some distribution in some model M and P is some structure inducing
convex function, typically of the form (10) or (11). Cone constraints are also
easily handled.
The subgradient equation for such a problem at a solution β̂(S) reads
0 ∈ α̂(S) + ẑ(S)

(13)

(S, β̂(S), α̂(S), ẑ(S)) ∈ S F (ℓ, P),

(14)

with
where

def

S F (ℓ, P) =


(s, β, α, z) :

s ∈ supp(F ),
ℓ(β; s) < ∞,
α ∈ ∂ℓ(β; s),
P(β) < ∞,

z ∈ ∂P(β) .

(15)

Above, and throughout, ∂ and ∇ will denote subdifferentials and derivatives
with respect to β unless otherwise noted.
The set S F (ℓ, P) can be described by the “base space” supp(F ) and “fibers”
{(β, α, z) : ℓ(β; s) < ∞, α ∈ ∂ℓ(β; s), P(β) < ∞, z ∈ ∂P(β)} .
We call the map
θ̂

s 7→ (s, β̂(s), α̂(s), ẑ(s)) ∈ S F (ℓ, P)

(16)

the optimization map. In practice, given data S, a computer solves the problem,
i.e. produces a point in θ̂(S) ⊂ S F (ℓ, P). Specifically, a solver produces a point
in the range of the optimization map:

S0F (ℓ, P) = (s, β, α, z) ∈ S F (ℓ, P) : α + z = 0 .
(17)

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

11

Formally speaking, the program (12) may have no solutions. On this set θ̂(s) = ∅.
In all of our examples except the dual problem considered in Section B.4, we
will assume enough so that our convex programs have unique solutions, when
they have any.
The selection events Lee et al. (2016); Fithian et al. (2014) we are most interested in are typically of the form

SBF (ℓ, P) = (s, β, α, z) ∈ S0F (ℓ, P) : (β, α, z) ∈ B(s)
(18)

for some nice set-valued function B(s) which could be specified by the zero-set
of a function hB :
(β, α, z) ∈ B(s) ⇐⇒ hB (s, β, α, z) = 0.

In all examples below except basis pursuit in Section B.2, B(s) = B does not
depend on s. This set also has the form of a bundle with base space supp(F )
and fibers
{(β, α, z) ∈ B(s) : α + z = 0} .

We call SBF (ℓ, P) parametrizable if there exists a measurable parametrization
ψ defined on some domain D with range SBF (ℓ, P). The map ψ is typically
constructed to be an inverse of the optimization map θ̂ on SBF (ℓ, P).
Given a parameterization, we will typically
 construct a change of measure
using φ = πB ◦ ψ to simplify sampling from s : πB−1 (s) 6= ∅ by sampling from
SBF (ℓ, P) itself, where πB is the projection onto the base of SBF (ℓ, P). Our construction is similar to what Zhou (2014) called estimator augmentation. By
construction, then, πB (s, β, α, z) is such that
θ̂ ◦ πB = id|SBF (ℓ,P)
with id the identity map.
Transforming the probability space and constructing new data vectors that
solve (12) as functions of optimization problems is key to our approach of sampling. The explicit parametrization changes with each problem. Typically, in
selective inference examples, we condition on some function of β̂, perhaps its
support E and possibly the signs zE of the non-zero coefficients. More generally,
might condition on something besides (E, zE ) which we might denote by q. The
selection event, i.e. the quantity we condition on determines the set of constraints
Bq . We shall denote the parametrization map in our problems as ψq , so often we
will write ψ(E,zE ) . The map ψq always produces a point (s, β, α, z) ∈ SBFq (ℓ, P),
while the domain of this parametrization can vary with the problem at hand,
that is
ψq : Dq → SBFq (ℓ, P),
for a problem specific domain Dq .

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

12

2.1. Inverting the MLE
We now begin to describe our approach to inverting the optimization map. Our
first example is classical: the density of the MLE of the natural parameters of
an exponential family. The formula is not new, going back at least to Fisher
(Barndorff-Nielsen, 1983; Efron and Hinkley, 1978). Nevertheless, it serves to
illustrate the general approach we take for the general cases later.
Our loss function is
ℓ(β; S) = Λ(β) − β T S

where

eΛ(β) = EF0 [eβ

T

S

]

is the moment generating function in the exponential family with reference
measure F0 and sufficient statistic S with Fβ the law of S above and in what
follows. We also assume that F0 has a density f0 with respect to Lebesgue
measure.
As we are computing the MLE, our penalty function is P(β) ≡ 0. The convex
program we solve is
minimize
Λ(β) − β T S.
(19)
p
β∈R

The KKT conditions or subgradient equation here is just the usual score equation
∇Λ(β̂(S)) = S

where

∇Λ(β) =

Z

sfβ (ds).

Rp

Our parameterization of S0F (ℓ, P = 0) is
ψ(β) = (∇Λ(β), β, 0, 0)
which can be interpreted as reconstructing S given β.
Standard multivariate calculus then tells us that the density of β̂(S) under
the distribution Fβ0 is
e∇Λ(β)

T

β0 −Λ(β0 )

det(∇2 Λ(β))f0 (∇Λ(β)).

(20)

where ∇2 Λ(β) is the observed information (Efron and Hinkley, 1978).
If β0 is the true parameter, then, at the cost of changing the reference measure
by a factor of exp(sT β̂(s) − Λ(β̂(s)) det(∇2 Λ(β))1/2 yielding a new Lebesgue
density h0 , we can rewrite this as
eΛ(β)−Λ(β0 )+∇Λ(β)

T

(β0 −β)

det(∇2 Λ(β))1/2 h0 (∇Λ(β)).

(21)

As pointed out in (Barndorff-Nielsen, 1983) ignoring the term h0 (∇Λ(β)) (which
is the constant 1 in the Gaussian case) yields the usual saddle-point approximation to the density of the MLE, up to the constant of integration. The exponential above can be rewritten as
1
− (β − β0 )T ∇2 Λ(β)(β − β0 ) + R(β; β0 ).
2
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

13

Note that even ignoring the remainder, this is quadratic in β0 the parameter,
and not β the variable of integration in the density.
In principle, nothing above really relies on the exponential family structure
for the model, though it does rely somewhat on the fact that the loss we use
came from an exponential family. It relies on this in that we use the form of
the loss to reconstruct data S from optimization variables β. This is similar to
what we see in the LASSO example below.
Nevertheless, the same argument above shows that if we solve the program
(19) then, so long as S ∼ F has a Lebesgue density the density of β̂(S) is
f (∇Λ(β)) · det(∇2 Λ(β)) .

(22)

In this sense, the above display provides an exact recipe to compute the density
of the MLE under model misspecification. This is somewhat similar to the general approach taken in Hillier and Armstrong (1999), though we are considering
this only in a very restricted setting.
2.2. Pull-Back of the LASSO with fixed design matrix
As a second example of a pull-back, we look at the canonical example in the class
of regularized convex optimization problems: the LASSO Tibshirani (1996b).
The LASSO program is defined for each (X, y, λ) ∈ Rn×p × Rn × (0, ∞) as
1
ky − Xβk22 + λkβk1 .
minimize
β∈Rp
2

(23)

In this example, X is considered fixed and F = L(PC y|X) where PC is projection onto col(X). The law F is supported on col(X) because the optimization
map depends only on PC y. In our general notation, we can take
S = PC y
1
ℓ(β; y) = ky − Xβk22
2
P(β) = λkβk1
n
S F (ℓ, P) = (y, β, α, z) : PC y = y, β ∈ Rp ,

(24)

o
α = X T Xβ − X T y, z ∈ ∂(λk · k1 )(β) .

The familiar subgradient equations of the LASSO are
− α̂(y) = X T (y − X β̂(y)) = ẑ(y),

ẑ(y) ∈ ∂(λk · k1 )(β̂(y)).

(25)

Note also that (25) contains the implicit constraint ẑ(y) ∈ row(X) where row(X)
is the rowspace of X. This can be seen from the structure of


F
F
T
S0 (ℓ, P) = (y, β, α, z) ∈ S (ℓ, P) : X (Xβ − y) + z = 0
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

14

as each point in S0F (ℓ, P) has z ∈ row(X).
In Lee et al. (2016), the authors assume general position so the map θ̂ is
single valued Tibshirani (2013). The authors then considered the active set and
signs of the LASSO solution, conditioning on their value (E, zE ) and the design
X. This is the event
n
o
y : (y, β̂(y), α̂(y), ẑ(y)) ∈ S0F (ℓ, P), diag(zE )β̂E (y) > 0, β̂−E (y) = 0 . (26)
Note that this event is equivalent to

θ̂(y) ∈ (y, β, α, z) ∈ S0F (ℓ, P) : diag(zE )βE > 0, β−E = 0 .
This is our canonical example of a parameterizable set where

B = B(E,zE ) = {(β, α, z) : diag(zE )βE > 0, β−E = 0, zE = λ sign(βE ), kz−E k∞ ≤ λ} .
The above constraints on β induce further restrictions on the range of the solver
SBF(E,z ) (ℓ, P) = {(s, β, α, z) ∈ S0F (ℓ, P) : β ∈ B(E,zE ) }
E

= {(s, β, α, z) : diag(zE )βE > 0, β−E = 0,

α = −X T (s − Xβ) = −z, zE = λ sign(βE ), kz−E k∞ < λ}.

The authors in Lee et al. (2016) then carry out selective inference for linear
functionals η T µ in the saturated model

(27)
M = N (µ, σ 2 ) : µ ∈ Rn

with σ 2 > 0 considered known. In this context, selective inference corresponds
to taking each F ∈ M and conditioning it on the event (26) which can be
rewritten as


πS SBF(E,z ) (ℓ, P)
E

where πS denotes projection onto the data coordinate. The resulting selective
model

)
(
(
F
∗
S
(ℓ,
P)
1
y
∈
π
dF
S
B
(E,zE )
(y) ∝
M∗ = F ∗ :
dF
0 otherwise.

is an exponential family and sampling is generally not necessary in the saturated
model as valid inference typically requires conditioning on sufficient statistic
related to nuisance parameters as described in Lee et al. (2016); Fithian et al.
(2014).
Nevertheless, if sampling were necessary, a naive accept reject sampling scheme
for inference in this setting draws vectors y on Rn according to N (µ, σ 2 ) solves
the LASSO with the triple (y, X, λ) and checks whether the result lies in (26).
Due to the nature of the LASSO, this check can be reduced to verifying whether
y satisfies a set of affine inequalities Lee et al. (2016).

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

15

2.2.1. Parametrization and pull-back of the LASSO
What if we did not have to check these affine inequalities in our sampling
scheme? This is the essence of what we propose in this work. We will ultimately
see that our approach is essentially equivalent to that of Lee et al. (2016) but
the probability space of our sampler is different. The approach we take is similar to (Zhou, 2014). We will see later that, after randomization, the parameter
space is generically simpler than if we had not randomized.
Our first example of a parameterization is
ψ(E,zE ) (β, z) = (Xβ + (X T )† z, β, −z, z)

(28)

with domain
D(E,zE ) = {(β, z) : diag(zE )βE > 0, β−E = 0,

z ∈ row(X), zE = λ sign(βE ), kz−E k∞ ≤ λ}

(29)

and corresponding data reconstruction map φ(E,zE ) : D(E,zE ) → col(X) defined
by
φ(E,zE ) (β, z) = Xβ + (X T )† z.
We see that ψ(E,zE ) is indeed a parameterization of S0F (ℓ, P). Hence, on
the range of φ(E,zE ) we do not actually have to check the affine inequalities of
(Lee et al., 2016) as they are satisfied by construction.
The set D(E,zE ) is a subset of a 2p-dimensional set but is in fact of dimension
rank(X) whenever |E| < rank(X). Hence, it has Lebesgue measure 0 in R2p but
it is a subset of an affine space so that it inherits a Hausdorff measure Hrank(X) .
In a formal sense, the map φ(E,zE ) is the inverse of the optimization map.
Lemma 1. On the set of y where that the solution to (23) is unique (X being
considered fixed), the following equality holds
θ̂−(y,α) ◦ φ(E,zE ) (β, z) = (β, z) for (β, z) ∈ D(E,zE ) ,

(30)

where θ̂−(y,α) (y) = (β̂(y), ẑ(y)). Hence, θ̂−(y,α) ◦ φ(E,zE ) is equivalent to the
identity map on D(E,zE ) .
More generally, the set (26) is equal to φ(E,zE ) (D(E,zE ) ).
Proof. Direct verification of (25).
Now choose an element in M, i.e. fix some µ ∈ Rn and consider the N (µ, σ 2 I)
density. We define its pull-back to be the law Fµ∗ supported on D(E,zE ) with
Hrank(X) density


1
(31)
fµ∗ (β, z) ∝ exp − 2 kφ(E,zE ) (β, z) − µk22 · Jφ(E,zE ) (β, z) ,
2σ
where Jφ(E,zE ) is the Jacobian of the projection of the parameterization, which
depends only on X in this case (and can be treated as a constant). We use the
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

16

name pull-back as this density is precisely the pull-back of the measure with
Lebesgue density fµ under the reconstruction map φ(E,zE ) .
As each density has a pull-back, the model itself has a pull-back M∗ comprising the pull-back of each density. The model also inherits the exponential
family structure of M. Taking f0∗ to be the reference measure, we can choose
the natural parameter to be µ/σ 2 and sufficient statistic to be φ(E,zE ) (β, z).
Finally, note that nothing about our construction so far has anything to do
with the parametric model (27). In fact, as long as F has a Lebesgue density
the pull-back is well-defined.
Lemma 2. Suppose F = L(PC y|X) has density f on C = col(X) and the
solution to (23) is F -a.s. unique. Then, the following equality holds


Ly∼F PC y (y, β̂(y), α̂(y), ẑ(y)) ∈ S0F (ℓ, P), diag(zE )β̂E (y) > 0, β̂−E (y) = 0, X

= L(β,z)∼F ∗ φ(E,zE ) (β, z)|X ,
(32)
where F ∗ is supported on relint(D(E,zE ) ) and has Hrank(X) density

f ∗ (β, z) ∝ f φ(E,zE ) (β, z) · Jφ(E,zE ) (β, z) .
Proof. Follows from Lemma 1 and standard multivariate calculus combined with
the fact that the existence of the density f implies that
F (φ(E,zE ) (D(E,zE ) \ relint(D(E,zE ) ))) = 0.
Remark 1. One of the consequences of the above is that one can draw response
vectors y that have the same active set and signs from a density with relatively
simple support when row(X) = Rp . Hence, sampling IID from some density
on this support allows one to use importance sampling with the above explicit
density as numerator in the importance weight.
When row(X) ( Rp , the support is still somewhat complex. For instance, the
subgradients z must be in the intersection of a face of the ℓ∞ ball which may be
a difficult set to describe. We will see that after randomization, this complexity
often disappears.
Remark 2. In theory, one might want to drop the uniqueness assumption above.
If uniqueness of the solution does not hold then the parameterization is not
injective and the event we condition on should be replaced with the event θ̂(y) ∩
SBF (ℓ, P) 6= ∅.
In this case, it may still be possible to derive a formula for the law of PC y by
invoking the co-area formula Federer (1959); Diaconis et al. (2013). We do not
pursue this generalization here.
In order to carry out selective inference for a linear functional η T µ, the authors of Lee et al. (2016) conditioned on Pη⊥ y. This is also possible in the pullback model. Fix a linear subspace L ⊂ col(X) considered to be the model
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

17

subspace so that

ML = N (µ, σ 2 I) : µ ∈ L .

For each linear function of interest η ∈ L, one generally must condition on
(PL − Pη )y to eliminate the nuisance parameter (PL − Pη )µ. If w is the observed
value of (PL −Pη )y, a straightforward modification of the argument above shows
that the appropriate distribution is supported on the rank(X) − dim(L) + 1
dimensional set
D(E,zE ) (w) = {(β, z) ∈ D(E,zE ) : (PL − Pη )φ(E,zE ) (β, z) = w}
with Hrank(X)−dim(L)+1 density proportional to fµ∗ . In this case, even if row(X) =
Rp , the support is somewhat complex.
3. Inverting the optimization map of a randomized convex program
In this section, we describe how to invert the KKT conditions for our family of
randomized convex programs. As mentioned in the introduction, randomization
comes with advantages like enhanced statistical power in the inference stage,
as well as a simplification of the support of the relevant reference distribution.
In what follows, all of our convex programs have random variables appearing
linearly in the subgradient, so that these random variables can be reconstructed
from optimization variables. The random variable is one introduced by a data
analyst through additional randomization as described in (Dwork et al., 2015;
Tian and Taylor, 2015).
Let G be a distribution on Rp and ǫ ≥ 0 a small parameter. Given a canonical
problem specified by (F, ℓ, P) we define its randomized version as
ǫ
minimize
ℓ(β; S) + P(β) − ω T β + kβk22 ,
β∈Rp
2

(S, ω) ∼ F × G.

(33)

A natural question to ask at this point is: why randomize the program in
the above way? One of the inspirations for selective inference after randomization, the topic considered in (Tian and Taylor, 2015) are techniques used in
differential privacy (Dwork et al., 2015). The other was the empirical finding in
(Fithian et al., 2014) that holding out some data before carrying out selective
inference as in Lee et al. (2016) also showed an improvement in selective power.
This increase in selective power can be attributed to the fact that choosing
a model with a randomized response has more leftover information as defined
in Fithian et al. (2014) after model selection than choosing a model with the
original response.
The choice to add 2ǫ kβk22 − ω T β to the objective is not the only reasonable choice. The main property we want of the perturbation ∆(β, ω) is that
ω = φ(β, v) where (β, v) ∈ ∂∆(β, ω). That is, given the value of β and the
subgradient v we can reconstruct ω.

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

18

Remark 3. Addition of the term 2ǫ kβk22 to the objective ensures that the set of
(s, ω) such that (33) has a solution contains the set
{(s, ω) : θ̂(s) 6= ∅, ω ∈ supp(G)}.
To see this, note that solving (33) is equivalent to computing
β̂(s, ω) = prox 1ǫ (ℓ(·;s)+P(·)) (ω/ǫ).
Whenever θ̂(s) 6= ∅, the objective ℓ(·; s)+P(β) is a closed proper convex function,
hence its proximal mapping is well-defined and β̂(s, ω) is well-defined.
If ℓ is strongly convex a.s.-F and P is closed and proper then we can take
ǫ = 0.
3.1. Data splitting as randomization
The suggested randomization above may strike some readers as lacking motivation. A randomization mechanism probably familiar to most readers is data
splitting (Cox, 1975), i.e. randomly partitioning the dataset into two pieces of
size (n1 , n − n1 ). We will see that solving a convex program after data splitting
is essentially an example of (33).
For a specific example, suppose we fit the graphical LASSO Friedman et al.
(2008) based on a data matrix X ∈ Rn×p but we first randomly split the data.
That is, we form


1
1
T
T
S1 (g) =
X1 (g),
In1 −
X1 (g)
1n 1
n1
n1 1 n1
with g denoting the random partition of the data into two groups and X1 (g)
denoting the data in the first group.
Next, we solve the program
minimize − log det(Θ) + Tr(S1 (g)Θ) + λP(Θ)

Θ:ΘT =Θ,Θ>0

with the usual caveat that the penalty does not charge diagonal elements of Θ.
This program is equivalent to solving the problem
minimize − log det(Θ) + Tr(SΘ) + Tr((S1 (g) − S)Θ) + λP(Θ),

Θ:ΘT =Θ,Θ>0

where
S=

1 T
X
n



1
I − 1n 1Tn X.
n

This is very close to (33) with ω = ω(g) = S − S1 (g) which is orthogonal
to, but not necessarily independent of the full covariance S. As described in
Fithian et al. (2014), data splitting can be interpreted as simply conditioning
on the value (S1 (g), g) as the only variation in L(S|S1 (g), g) is S2 (g), the second
stage data.
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

19

3.2. The selective sampler
In this section, we describe our main tool for inference after solving a convex
program randomized in the above fashion. We call this tool the selective sampler.
A solver for the program (33) produces a point in

def
S̄ F (G, ǫ, ℓ, P) = (s, ω, β, α, z) :
s ∈ supp(F ),

ω ∈ supp(G),
ℓ(β; s) < ∞,

(34)

α ∈ ∂ℓ(β; s),
P(β) < ∞,

z ∈ ∂P(β) .

More precisely, it produces a point in

S̄0F (G, ǫ, ℓ, P) = (s, ω, β, α, z) ∈ S̄ F (G, ǫ, ℓ, P) : ǫ · β + α + z − ω = 0 . (35)

Our statistical learning task is typically to infer something about F in some
model M. As we are free to choose G however we want, a natural choice is
to choose G to have a Lebesgue density supported on all of Rp . In this case,
inspection of the KKT conditions or subgradient equation of (33) read
ω = α̂(s, ω) + ẑ(s, ω) + ǫ · β̂(s, ω).

(36)

The subgradient equation determines a canonical map
ψ : S F (ℓ, P) → S̄0F (G, ǫ, ℓ, P)

(37)

ψ(s, β, α, z) = (s, α + z + ǫ · β, β, α, z).

(38)

defined naturally as

As in the non-randomized case, we are typically interested in some selection
events of the form
S̄BF (G, ǫ, ℓ, P) = {(s, ω, β, α, z) ∈ S̄0F (G, ǫ, ℓ, P) : (β, α, z) ∈ B(s)}.
The corresponding randomization reconstruction map is then
φ(β, α, z) = α + z + ǫ · β,
defined on a domain
{(β, α, z) ∈ B(S) : α + z + ǫ · β ∈ supp(G)},
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

20

for a suitable B(S), defined by the selection event. It turns out that many
selection events of interest are such that the restriction of ψ to these events have
a simple structure which allows for straightforward description of the selective
model. The canonical example of conditioning on the set of active variables and
signs of the LASSO as in (Lee et al., 2016) was described in Section 1.1.
Theorem 1 (Selective sampler). Suppose that ω is independent of S with distribution G such that supp(G) ⊂ Rp has non-empty interior with Lebesgue density
g on supp(G). Then, the map
ψ(s, β, α, z) = (s, ǫ · β + α + z, β, α, z) = (s, φ(β, α, z), β, α, z)
restricted to

DG = (s, β, α, z) ∈ S F (ℓ, P) : ǫ · β + α + z ∈ supp(G)

is onto S̄0F (G, ǫ, ℓ, P). Further, the law

LF ×G ((S, ω)|(S, ω, β̂(S, ω), α̂(S, ω), ẑ(S, ω)) ∈ S̄BF (G, ǫ, ℓ, P))
= L((S, ǫβ + α + z)|(s, β, α, z) ∈ DG , (β, α, z) ∈ B(S))

(39)

for suitable B(S) and (S, β, α, z) has density proportional to
f (s) · g(ǫ · β + α + z) · |Jψ(s, β, α, z)| · 1DG (s, β, α, z) · 1B(s) (β, α, z)

= f (s) · g(φ(β, α, z)) · det(D(β,α,z) φ) · 1DG (s, β, α, z) · 1B(s) (β, α, z)

(40)

with the Jacobian denoting the derivative of the map ψ with respect to (β, α, z)
on the fiber over s.
Proof. The fact that ψ restricted to DG is onto follows from its construction. Let
F (s) = {(β, α, z) : α ∈ ∂ℓ(β; s), z ∈ ∂P(β)} denote the fiber over s and assume
that B(s) = B does not vary with s. In this case, in any local coordinates on
F (s) standard multivariate calculus yields the density (40) as the derivative of
the map
(s, β, α, z) 7→ (s, α + z + ǫ · β)
takes the form


In×n
Ds φ

0n×p



D(β,α,z) φ

with determinant det(D(β,α,z) φ). The result then follows by integrating over
supp(F ) with density f (s). If S does not have a density then (39) can be derived
via the Kac-Rice formula Taylor et al. (2013); Adler and Taylor (2007).
If B depends on s then the Kac-Rice formula may often be used to derive
the above density given the B can be described as the zero set of some smooth
function hB . As most of our examples are such that B does not depend on s we
omit the details.

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

21

Remark 4. The correct interpretation of (39) has as conditioning event
θ̄−(s,ω) (s, ω) ∩ B(s) 6= ∅
where θ̄ is the map
θ̄(s, ω) = (s, ω, β̂(s, ω), α̂(s, ω), ẑ(s, ω)).
Remark 5. When supp(G) = Rp then DG = S F (ℓ, P).
Remark 6. Above, we have constructed ω as a function of (s, β, α, z). This
is similar to the LASSO case in Section 2.2 in that we construct new random
variables out of optimization variables. Often, it is also possible to construct s
from (ω, β, α, z) though there is no canonical embedding unless we provide more
structure to the map ℓ. If ℓ is an exponential family negative log-likelihood so
that
ℓ(β; s) = Λ(β) − sT β
then this is certainly possible. However, in order to have a change of measure
result as in Theorem 1, the law of S should have a density on Rn .
Remark 7. Often, we will want to condition on some functions of S. Theorem 1
formally holds unchanged for any distribution supported on a lower dimensional
subset of S. One simply replaces the law F with the appropriate law supported on
a lower dimensional set. This device was used in Section 1.1. A further example
of this is considered in Appendix A below.
Remark 8. The theorem assumes ω is independent of the data S. It is straightforward to see that a similar result holds if we replace throughout the density g
above with K(ω; s) a kernel for the conditional density of ω|s.
4. Polyhedral examples
We now begin describing several instances of the selective sampler. In this section, the penalties or constraints are polyhedral. In this case, the reconstruction
maps are typically affine in the optimization variables but may be non-trivial in
the data. Some of these examples were also considered in Tian et al. (2016). We
repeat them here, in more explicit detail, as concrete examples of the selective
sampler.
4.1. LASSO with Gaussian errors and fixed design matrix
As is often the case, the LASSO serves as a canonical example. We denote the
parametrization, based on the active set and signs as ψ(E,zE ) on domain D(E,zE )
and the reconstruction map for randomization ω as φ(E,zE ) throughout. In all
the below examples, we are implicitly thinking of cases when supp(G) = Rp and
G has a Lebesgue density, the canonical example being N (0, σω2 Ip ).
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

22

This example was addressed in Section 1.1, though we present it here in the
general notation developed so far. The embedding in Theorem 1 plays the role
of ψ(E,zE ) in the parametric LASSO example. The randomized LASSO program
Tian and Taylor (2015) with randomization G is defined for each (y, X, λ, ω) ∈
Rn × Rn×p × (0, ∞) × Rp as
ǫ
1
ky − Xβk22 + kβk22 − ω T β + λkβk1 ,
minimize
p
β∈R
2
2

y|X ∼ F, ω ∼ G

(41)

with ω independent of (X, y).
When G = δ0 , we recover the parametric LASSO with ǫ = 0 and the parameteric Elastic Net if ǫ > 0 Zou and Hastie (2005). Supposing then that
supp(G) = Rp , and the canonical selection event given by
B(E,zE ) = {β : diag(zE )βE > 0, β−E = 0} .
A parametrization of
{(s, ω, β, α, z) ∈ S̄0F (G, ǫ, ℓ, P) : β ∈ B(E,zE ) }
is given by

 


βE
zE
T
ψ(E,zE ) (y, βE , u−E ) = y, ǫ
,
− X (y − XE βE ) + λ
0
u−E
 


βE
zE
, −X T (y − XE βE ), λ
0
u−E
with the domain
D(E,zE ) = {(y, βE , u−E ) : diag(zE )βE > 0, ku−E k∞ ≤ 1} .
The reconstruction map for ω is given by
 


β
zE
φ(E,zE ) (y, βE , u−E ) = ǫ E − X T (y − XE βE ) + λ
,
0
u−E

(42)

(43)

again with the same support as the parametrization map. For the canonical
event B(E,zE ) , we therefore need to sample from a density proportional to
  


β
zE
T
fµ (y) · g ǫ E − X T (y − XβE ) + λ
· det(XE
XE + ǫI)
(44)
0
u−E
supported on D(E,zE ) , where fµ is the N (µ, σ 2 In ) density. As above, the same
expression holds if our model for y|X is not from the normal family.
A common variant of the LASSO implemented in Friedman et al. (2010) is
minimize
p
β∈R

ky − Xβk22
+ λkβk1 .
2kX T yk∞

(45)

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

23

A randomized version of this is easily handled, changing the density to be proportional to


 T

  
XE XE
X T (y − XE βE )
zE
β
+
λ
·
det
+
ǫI
fµ (y) · g ǫ E −
u−E
0
kX T yk∞
kX T yk∞
and supported on the same set as in (42). Another variant Sabourin et al. (2014)
replaces kX T yk∞ with a resampling based quantity median(X T y) resulting in
a similar change in the sampling density.
Remark 9. Conditioning on Nuisance Statistics: As detailed in Fithian et al.
(2014), we condition not just on the selection event B(E,zE ) , but also on the sufficient statistics corresponding to nuisance parameters to obtain optimal UMPU
selective tests/ intervals. In such a case, the support for y in the selective sampling density is restricted to a set, denoted as Dobs . This has been sketched for
interested readers in A in the appendix.
4.2. The selective sampler is model agnostic: LASSO without
parametric assumptions
Up this point, we have assumed so far that X is fixed and the law of y|X was
from the parametric model {N (µ, σ 2 I) : µ ∈ Rn }.
We now remove this assumption, assuming that the law F is now just a law
for the pair (X, y) ∈ Rn×p × Rn . A common assumption is the pairs model:
IID
(xi , yi ) ∼ F̃ for some distribution F̃ on Rp × R, where xTi are the rows of X.
While we keep this model in mind, it is not necessary in order to define the
appropriate sampler. This section serves as an example of the sense in which
our main result Theorem 1 is agnostic to the underlying statistical model.
In the notation established so far
S = (X, y) ∈ Rn × Rn×p ,
1
ℓ(β; (X, y)) = ky − Xβk22 ,
2
P(β) = λkβk1 .
Therefore, we must sample from a density proportional to
  


β
zE
T
f (X, y) · g ǫ E − X T (y − XE βE ) + λ
· det(XE
XE + ǫI)
0
u−E
supported on
{(X, y, βE , u−E ) : diag(zE )βE > 0, ku−E k∞ ≤ 1} .

(46)

In order to use this result for particular applications of inference, one may
have to reduce the problem to sampling from a particular reference distribution perhaps by conditioning on appropriate sufficient statistics, as described in
Appendix A.
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

24

4.3. Selective CLT: ℓ1 -penalized logistic regression
In this section, we describe an application of the selective CLT of Tian and Taylor
(2015), applied to the case of logistic regression with random design matrix
X ∈ Rn×p with rows xTi , i = 1, . . . , n. Suppose
iid

xi ∼ FX , xi ∈ Rp , yi |xi ∼ Bernoulli(π(xTi b)),
where π(x) = exp(x)/(1 + exp(x)) and b ∈ Rp is unknown, p fixed, S = (X, y).
With logistic loss,
#
" n
1 X
T
T
yi log π(xi β) + (1 − yi ) log(1 − π(xi β)) ,
ℓ(β; (X, y)) = − √
n i=1
we solve randomized ℓ1 -penalized logistic regression introduced in Tian and Taylor
(2015),
ǫ
(47)
minimize
ℓ(β; (X, y)) − ω T β + kβk22 + λkβk1 .
β∈Rp
2
On the usual selection event of observing active set and signs (E, zE ), the randomization reconstruction map is


 
1 T
zE
βE
.
φ(E,zE ) (y, βE , u−E ) = ǫ
− √ X (y − π(XE βE )) + λ
u−E
0
n
And, the sampling density becomes proportional to

  

1
β
zE
f (X, y) · g ǫ E − √ X T (y − π(XE βE )) + λ
0
u−E
n


1
T
W (XE βE )XE + ǫI ,
· det √ XE
n

(48)

where W (Xβ) = diag(π(Xβ)(1 − π(Xβ)))1 is the Hessian matrix of the loss
and the density above is supported on the same set as in (46). Sampling (X, y)
jointly is not feasible when FX is unknown. Denoting
β̄E = argmin −
β∈R|E|

n
X
i=1

yi log π(xTE,i β) + (1 − yi ) log(1 − π(xTE,i β)),

(49)

the MLE for the unpenalized logistic regression with only the variables in E, a
Taylor series expansion of ∇ℓ(β; (X, y)) gives2


√
Q(βE − tE )
∇ℓ(β; (X, y)) ≈ n
.
−t−E + C(βE − tE )
1 With slight abuse of notation, we allow π : Rn → Rn , x 7→ π(x) to be the function applied
on each coordinate of x ∈ Rn .
2 Alternatively, we might take β̄ to be the one-step estimator in the selected model starting
E
from β̂E Taylor and Tibshirani (2016).

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

Here, t is the observed value of the random vector


β̄E
,
T = 1 T
n X−E (y − π(XE β̄E ))

25

(50)

1 T
1 T
XE W (XE β̄E )XE and C = X−E
W (XE β̄E )XE .
n
n
Since Q and C converge by the law of large numbers to fixed quantities, we can
either treat them as fixed or bootstrap them. For more details, see Tian and Taylor
(2015).
With p is fixed, pre-selected, T properly scaled is asymptotically normal



√
bE
D
→ N (0, Σ),
(51)
n T−
γ−E
Q=

where Σ is estimable from the data. When the selected parametric model is
correct (E ⊇ supp b) we note that γ−E = 0 and Σ is in fact block-diagonal
(Taylor and Tibshirani, 2016). When the selected parametric model is not correct, one can estimate this covariance nonparametrically, using bootstrap (see
Section 7). Since asymptotically T is from an exponential family with parameters bE , we could base inference on bE by sampling from the distribution of T
instead of (X, y).
In this case, the parametrization map becomes




  
√
Q(βE − tE )
zE
β
,
+λ
ψ(E,zE ) (t, βE , u−E ) = t, ǫ E + n
u−E
−t−E + C(βE − tE )
0
 

 

βE √
Q(βE − tE )
zE
, n
,λ
,
0
−t−E + C(βE − tE )
u−E
(52)
|E|
with the domain Rp × RzE × [−1, 1]p−|E| . Hence, the density we sample from is
proportional to




  
√
Q(βE − tE )
zE
β
,
(53)
+λ
f (t) · g ǫ E + n
u−E
−t−E + C(βE − tE )
0
  
bE
restricted on the domain of ψ(E,zE ) , where f here is the density of N
,Σ .
0
We note that, in order to construct a valid test about some linear functional of
bE we can always condition on the observed u−E and / or βE if we desire. The
upside to such conditioning is that the sampling problem becomes easier, with
the downside usually being a loss in selective power.
4.4. First step of forward stepwise
The Kac-Rice tests, described in Taylor et al. (2013) are based on the solution
of
maximize η T X T y, where y ∼ F,
(54)
η∈K

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

26

X ∈ Rn×p is a fixed design, K a convex set that can be stratified into smooth
disjoint manifolds and the process η T X T y is Morse for almost every y ∈ Rn . The
simplest example of such tests is one step of forward stepwise model selection,
in which case K is the ℓ1 ball of radius 1. In this case X will be usually centered
and scaled so that X T y corresponds to the marginal Z statistics for p different
simple linear regressions.
Inference in this broad class of problems, for the global null, is based on the
test statistic η ∗T X T y, where
η ∗ = argmax η T X T y,

(55)

η∈K

with K being the polar set of convex set C. Since, the null distribution of the
above test statistic is intractable, we could instead provide inference based on
the conditional law


L η ∗T X T y η ∗ = argmax η T X T y .
η∈K

With this brief description of inference based on Kac-Rice tests, we focus back
on the randomized versions of the Kac-Rice objective in (54). The randomized
objective is given by
maximize
η T (X T y + ω) − IK (η), where y × ω ∼ F × G,
p
η∈R

(56)

with ω ∈ Rp and penalty manifests as the characteristic function of convex set
K, that is
(
0 if η ∈ K,
IK (η) =
∞ otherwise.
Here we set ǫ = 0, as the above optimization problem does always have a
solution. Having ǫ > 0 would allow for several variables to be selected.
As mentioned above, perhaps the simplest example of the randomized optimization problem in (56) is forward stepwise problem with a single step. Specifically, that is consider
maximize η T (X T y + ω),
η∈K

for
K = {η ∈ Rp : kηk1 ≤ 1}.
The above optimization problem yields the optimal direction
(
s∗ if j = j ∗
∗
ηj =
0
otherwise,
j = 1, . . . , p, where
j ∗ = argmax |XjT y + wj |,
1≤j≤p

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

27

coordinate with the maximum absolute value, and
s∗ = sign(XjT∗ y + wj ∗ ),
the corresponding sign. Thus, we condition on the first active direction (both
s∗ and j ∗ ), which gives rise to selection event

Ê(s∗ ,j ∗ ) = (y, ω) ∈ Rn × Rp : sign(XjT∗ y + ωj ∗ ) = s∗ ,
s∗ (XjT∗ y + ωj ∗ ) ≥ max |XjT y + ωj | .
1≤j≤p

The subgradient equation yields the reconstruction map
φ(j ∗ ,s∗ ) (y, z) = z − X T y,
where z ∈ ∂IK (η ∗ ), the set of sub-gradients to K at η ∗ , given by the normal
cone
∂IK (η ∗ ) = {c(u1 , . . . , uj ∗ −1 , s∗ , uj ∗ +1 , . . . , up ) : ui ∈ R, |ui | ≤ 1, c > 0}.
We can reparametrize the set
n
o
(y, ω, β, α, z) ∈ S̄0F (G, ǫ = 0, ℓ, P) : (y, ω) ∈ Ê(s∗ ,j ∗ )

using

ψ(s∗ ,j ∗ ) (y, z) = (y, z − X T y, η ∗ , X T y, z),

(ηj∗ = s∗ I{j=j ∗ } , j = 1, . . . , p), with the domain Rn × ∂IK (η ∗ ).
With the above reparametrization, we sample (y, z) from a density proportional to
f (y) · g(z − X T y),
(57)
supported on Rn × ∂IK (η ∗ ).

Remark 10. The set of sub-gradients at η ∗ for s∗ > 0 can be identified as the
epigraph of the ℓ∞ norm (modulo a permutation of the maximum coordinate j ∗ )
and for s∗ < 0, ∂IK (η ∗ ) is the polar cone of the epigraph of ℓ1 norm.
Remark 11. Of course, it is usually of interest to take more than one step of
forward stepwise. Inference after several steps of forward stepwise is considered
in Tibshirani et al. (2014); Fithian et al. (2015). We consider several steps of
forward stepwise in Section 6.
Remark 12. The data analyst may have some set of variables Ē that she
insists on controlling for. In this case, instead of just assuming X is centered,
⊥
we might assume that PĒ
X = X and its columns normalized, where PĒ denotes
projection onto the column space of XĒ . In this way, centering X corresponds
to the common practice controlling for an intercept in the model.

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

28

Remark 13. Forward stepwise is also used in generalized linear models such
⊥
as logistic regression. In this case, the quantity X T PE
′ y that appears in the
least squares setting can be replaced either by the Wald Z-statistics or the score
Z-statistics. For example, for logistic regression, the score Z-statistics take the
form
XjT (y − π(XE ′ β̄E ′ ))
,
1≤j≤p
(58)
Zj =
(XjT W (XE ′ β̄E ′ )Xj )1/2
where β̄E ′ is the unpenalized MLE for the model with variables E ′ and
W (XE ′ β̄E ′ ) = diag(π(XE ′ β̄E ′ )(1 − π(XE ′ β̄E ′ ))
is a consistent estimate of the variance of y under the model with variables E ′ .
While this is the typical parametric estimate of variance forward stepwise would
use, one might prefer using a jackknife or bootstrap estimate of this variance if
one is unsure whether the model with variables E ′ is correctly specified, as would
be the case early on in building a model via forward stepwise. Using Wald type
Z statistics would require fitting p different logistic regression models which has
some computational burden.
Having computed the Zj ’s one might then consider a randomized version of
the problem
maximize η T Z
η:kηk1 ≤1

and proceed as above with S being the vector Z.
4.5. Marginal screening
Marginal screening computes marginal Z statistics
Zj =

XjT y
σkXj k2

Tj =

XjT y
σ̂j kXj k2

or T statistics

for each of p centered variables Xj and thresholds their absolute value at some
threshold, perhaps z1−α/2 where α is some nominal p-value threshold. This can
be expressed in optimization form as
1
minimize kη − T k22 .
η:kηk∞ ≤c 2
Selective inference in the nonrandomized setting for this problem was considered
in (Lee and Taylor, 2014).
A natural randomized version would be
1
minimize kη − T k22 − ω T η.
η:kηk∞ ≤c 2
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

29

Conditioning on the set achieving the threshold c and their signs to be (E, zE ),
we see that this event is
{(T, η, z) : ηE = c · zE , diag(zE )zE ≥ 0, kη−E k∞ < c, z−E = 0}

(59)

and the randomization reconstruction map becomes

 

z
c · zE
−T + E .
φ(E,zE ) (T, η−E , zE ) =
0
η−E
We thus sample (y, X, η−E , zE ) (or (T, η−E , zE ) if X is random) from a selective
density proportional to

 

z
c · zE
−T + E
f (T ) · g
(60)
0
η−E
and supported on the event in (59), where f is the unselective law of T . For
logistic regression, one can replace the T statistics above with the score statistics
as described in Remark 13.
4.6. Full model screening
Another possible way to screen variables is to threshold the coefficients from
the output of an unpenalized, randomized selection program given by
ǫ
1
ky − Xβk22 + kβk22 − ω T β, ((X, y), ω) ∼ F × G,
minimize
p
β∈R
2
2

(61)

Based on the solution β̂((X, y), ω), we define the threshold model (E, zE ) as
E = {i : |β̂i ((X, y), ω)| ≥ aσ}
and signs zE = sign(β̂E ((X, y), ω)), where a is a constant (e.g. 1 − α quantile
of the standard normal or the t-distribution for some nominal level α) and σ is
the scaling. The canonical event of interest is
B(E,zE ) = {β ∈ Rp : |βi | ≥ aσ ∀i ∈ E, kβ−E k∞ < aσ, diag(zE )βE > 0}.
A parametrization of
{((X, y), ω, β, α, z) ∈ S̄0F (G, ǫ, ℓ, P) : β ∈ B(E,zE ) }
is given by
ψ(E,zE ) ((X, y), β) = (X, y), ǫβ − X T (y − Xβ), β, −X T (y − Xβ), 0
with the domain



D(E,zE ) = {((X, y), β) : diag(zE )βE > 0, |βi | ≥ aσ ∀i ∈ E, kβ−E k∞ < aσ} .
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

30

For the canonical event B(E,zE ) , we therefore need to sample ((X, y), β) from
a density proportional to

f (X, y) · g ǫβ − X T (y − Xβ) · det(X T X + ǫI)
(62)

supported on D(E,zE ) , where f and g are the densities of F and G, respectively. If
X is random, the sampling of the data simplifies as in Section 4.3. The scaling σ
can be estimated as the variance of noise from the selected model with response
y and matrix of predictors XE . As long as the estimate of σ is consistent, we
can treat it as a constant in the selection event, hence have the selection event
as the polyhedral region.
4.7. Selective sampler based on a penalized MLE: Fisher’s exact
selective test
In Section 2.1, we considered the pull-back of the unpenalized MLE for an
exponential family. Allowing for penalization in the problem (19) as well as
randomization yields programs of the form
ǫ
minimize
Λ(β) − S T β + kβk22 − ω T β + P(β), (S, ω) ∼ F × G.
β∈Rp
2
We describe a version of Fisher’s exact test based on observing some function
of the solution to the above program. As we saw in the LASSO, the penalty P
is often chosen such that the solution set
{(β, z) : z ∈ ∂P(β)}
possesses a nice stratification into a discrete collection of bundles C, the canonical
example being the ℓ1 norm, or perhaps non-negative constraints. In particular,
if P is the support function of some convex set K, then (4.7) is seen to be the
normal bundle of K which itself often has a nice stratification.
Fixing C to be one of these strata, the subgradient equations of the optimization problem yields a reconstruction map
φC (s, β, z) = ∇Λ(β) − s + z + ǫ · β
where (β, z) ∈ C and s ∈ supp(F ). In this case, the parametrization map of
{(s, ω, β, α, z) ∈ S̄0 (G, ǫ, ℓ, P) : (β, z) ∈ C}
becomes
ψC (s, β, z) = (s, ∇Λ(β) − s + ǫβ + z, β, ∇Λ(β) − s, z)
and the sampling density proportional to
f (s) · g(∇Λ(β) − s + ǫβ + z) · |JφC (s, β, z)|
and restricted to (s, β, z) ∈ M × C.
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

31

We now describe a stylized instance of such a problem used in the context of
goodness-of-fit tests. Much of the effort in selective inference has been focused
on finding tools for inference about specific parameters in a model, i.e. inference
about statistical functionals θ : M → R. Two examples in the literature that
employ such goodness-of-fit tests are Choi et al. (2014); Fithian et al. (2015).
In this case, we assume the data analyst will use the data to decide what
sufficient statistics to use in constructing the goodness-of-fit test. Suppose our
data S is a poisson process N with points in X and intensity measure modeled
with density


p
X
dλ(β)
βj hj (x)
(x) = exp 
dµ
j=1

with respect to some reference probability measure µ on X .
If X = [0, 1] and hj = 1[lj ,uj ] were a collection of indicator functions, this
would correspond to a change-point model. We might then use something like
the fused LASSO in which we can take the intervals to be the step functions
[lj , uj ] = [(j − 1)/p, 1] or perhaps a multiscale approach as in Frick et al. (2014);
Chan and Walther (2011); Rivera and Walther (2013) and use the penalty
P(β) =

p
X
j=1

λj |βj |

a weighted ℓ1 penalty3 . Usually, the background rate is included in such a model
by fixing h1 to be the constant function 1. An unpenalized fit for this variable
corresponds to its corresponding λ being 0.
Depending on our choice of penalty, we might then solve the following randomized program

Z
p 
X
ǫ
β
minimize
Λ(β)
−
h
(x)
N
(dx)
+ P(β) + kβk22 − ω T β
j
j
p
β∈R
2
X
j=1
where



Λ(β) = log 

Z

X





exp 

p
X
j=1







βj hj (x) − 1 µ(dx) .

Having solved the above problem, the data analyst now observes that (β̂(N, ω), ẑ(N, ω))
are in some subbundle C of the set of variable-subgradient pairs. When P is a
weighted ℓ1 penalty, the conditioning event
n
o
(N, ω) : Ê(N, ω) = E, zÊ(N,ω) = zE
where

n
o
Ê(N, ω) = j : β̂j (N, ω) 6= 0, λj > 0

3 The multiscale approach described in Frick et al. (2014); Chan and Walther (2011);
Rivera and Walther (2013) is formally a testing approach. The penalized version above was
proposed in Tian et al. (2015) based on the dual of the test statistic in the multiscale literature.

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

32

are the selected penalized coordinates. The corresponding subbundle is the set

C = (βE , βU , z−E ) : diag(zE )βE ≥ 0, kz−E k∞ ≤ 1, βU ∈ RU

where U = {j : λj = 0} is the set of unpenalized coordinates. With some abuse
of notation, we write z = (zE , zU , z−E ) = (zE , 0, z−E ) as the full subgradient
where zE are held fixed on C and zU = 0 as these coordinates are unpenalized.
The data analyst now decides to test the null hypothesis
H0,j|E : βj|E = 0,
in the model
log λ(β) =

X

βj

Z

j 6= 1
hj (x) µ(dx).

X

j∈E

Following the exponential family setup in Fithian et al. (2014), we might condition on the nuisance sufficient statistics
Z
hi (x)N (dx),
i ∈ E \ j.
X

Assuming we have included a background rate in the model, this fixes the total
number of points in the Poisson process to be N (X ) the observed number of
points.
If we had not used the data to select the intervals, conditioning on these
sufficient statistics and resampling points is exactly Fisher’s exact test modulo
the choice of test statistic. That is, the appropriate reference measure can be
constructed by sampling from a Binomial process Ñ with n = N (X ) points and
distribution µ conditioned to satisfy
Z
Z
hi (x) N (dx),
i∈E\j
hi (x) Ñ (dx) =
X

X

where the right hand side are the values observed in the data. When X is
a discrete space, this is a generalization of Fisher’s exact test. Sampling for
such problems have a fairly rich literature (c.f. (Diaconis and Sturmfels, 1998)).
When X is Euclidean or a manifold, the conditioning event above corresponds to
a subset of the configurations of N (X ) points on X and is generally a nontrivial
task. See (Diaconis et al., 2013) for further discussion and examples.
Having used the data to choose which sufficient statistics to use, we must
use an appropriate selective distribution. Beyond just sampling Ñ from the
conditional density, we must sample (βE , βU , z−E ). Conditional on Ñ we see
that the appropriate density of the joint law with respect to the product of the
law of Ñ and p-dimensional Hausdorff measure on C is proportional to


Z
(Ñ , βE , βU , z) 7→ g ∇Λ(β) −
h(x)dÑ (x) + ǫ · β + diag(λ)z
X

with
∇Λ(β) =

R


P
βi hi (x) µ(dx)

P i∈E∪U
exp
i∈E∪U βi hi (x) µ(dx)

XRh(x) exp
X

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

33

For computational reasons, to evaluate the integrals over X above, it may be
simpler to use a discretization of X as in Lindsey’s method Efron and Tibshirani
(1996).
Finally, while we have considered using a weighted LASSO to choose the sufficient statistics, one might use a penalty with some curvature as well, requiring
the modifications discussed in Section 5.
4.8. Graphical models
Gaussian graphical models are a popular way to study network structures. In
particular, it has often been used on many types of genome data, e.g. gene expression, metabolite concentrations, etc. Specifically, consider the p-dimensional
normally distributed random variable
X = (x1 , . . . , xp ) ∼ N (µ, Σ).
It is of interest to study the conditional independence structure of the variables
{1, 2, . . . , p}. The conditional independence structure is conveniently represented
by an undirectional graph (Γ, E), where the nodes are Γ = {1, 2, . . . , p}, and
there is an edge (i, j) ∈ E if and only if xi 6⊥ xj conditional on all the other
variables Γ\{i, j}. Moreover, assuming the covariance matrix Σ is not singular,
we denote the inverse covariance matrix Θ = Σ−1 , then
xi ⊥ xj |XΓ\{i,j} ⇐⇒ Θij = 0.
In many applications of Gaussian graphical models, we assume the sparse edge
structure, where we can hope to uncover the network structure even in the highdimensional setting. We discuss applying the selective sampler to the graphical
LASSO Friedman et al. (2008). Another graphical model selection, neighborhood selection, is discussed in Appendix C.
The randomized graphical LASSO with randomization ω ∼ G is the convex
program
minimize − log det(Θ) + Tr((S − ω)Θ) + λP(Θ),

Θ:ΘT =Θ,Θ>0

where
P(Θ) =

X

(i,j):i6=j

(S, ω) ∼ F × G, (63)

|Θ|ij

and
Usually,


S ∈ Mk = A : AT = A, A ≥ 0, rank(A) = k ,
S=

1 T
X
n



1
I − 11T X,
n

k ≤ p.

X ∈ Rn×p

for some data matrix so that k is generically min(n − 1, p). Note that Θ, S, ω
are p × p matrices.
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

34

The subgradient equations read as
ω = Θ̂−1 (S, ω) + S + ẑ(S, ω),
where ẑ(S, ω) ∈ ∂P(Θ̂(S, ω)). For E ⊂ {(i, j) : i 6= j}, the natural selection
event of interest is
n
o
S = (S, ω) : ∃ Θ = Θ̂(S, ω) s.t. sign(ΘE ) = zE , Θ−(E∪{(i,i):1≤i≤p}) = 0 ,

Owing to symmetry in Θ, on S we have

2λsij
∂P(Θ̂(S, ω)) 
ẑij (S, ω) =
= 0

∂Θij

2λuij

(i, j) ∈ E
i=j
(i, j) ∈
/ E and i 6= j,

where
and

uij = uji satisfies kuk∞ ≤ 1 for (i, j) ∈
/ E, i 6= j,
sij = sign(Θij ).

Due to symmetry, we can parametrize ω by the active elements in Θ along with
its diagonal elements (dimension |E|/2 + p) and the inactive elements of the
penalty subgradient (dimension (p2 − p − |E|)/2), corresponding to the inactive
elements, both restricted to the upper triangular part of Θ. To fix notations,
denote
E ∗ = {(i, j) ∈ E : i ≤ j} ∪ {(i, j) : i = j},
the set of active coordinates in the upper triangular part of Θ along with its
diagonal elements.
This leads us to the following parameterization of {(S, ω, Θ, α, z) ∈ S̄0 (G, ǫ =
0, ℓ, P) : (S, ω) ∈ S}
ψ(E,zE ) (S, ΘE ∗ , u−E ∗ ) = (S, Θ−1 (E ∗ ) + S + z(u−E ∗ ), Θ(E ∗ ),
Θ−1 (E ∗ ) + S, z(u−E ∗ )),

(64)

where Θ(E ∗ ) is the matrix with E ∗ entries from ΘE ∗ and zeroes elsewhere,
zij (u) = 2λ(sij I{(i,j)∈E} + ui,j I{(i,j)6∈E,i6=j} ) and the domain of φ(E,zE ) is restricted to ku−E ∗ k∞ ≤ 1.
The randomization reconstruction map can be written as
φ(E,zE ) (S, ΘE ∗ , u−E ∗ ) = Θ−1 (E ∗ ) + S + z(u−E ∗ ).
With this parametrization, we can sample from the density proportional to
f (S) · g(Θ−1 (E ∗ ) + S + z(u−E ∗ )) · |Jψ(S, ΘE ∗ , u−E ∗ )| .
and restricted to ku−E ∗ k∞ ≤ 1. The computation of above Jacobian involves
(see the Appendix of Taylor and Tibshirani (2016))
!
∂φ(E,zE )
= 2(σil σjk + σik σjl ) for (i, j) ∈ E ∗ ,
∂Θij
kl

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

∂φ(E,zE )
∂uij

!

kl

35

= 2λδik δlj for (i, j) ∈
/ E ∗ , i ≤ j,

where σij are the elements of Σ and δij is Kronecker delta function.
The above computations give the (k, l) coordinate of the (i, j)th upper triangular matrix and as again, we simply restrict attention to {(k, l) : k ≤ l}, the
upper triangular half of randomization map.
More polyhedral examples, as well as the generalized LASSO are presented
in Appendix B.
5. Curvature
Up to now, our examples have been polyhedral in nature. In most examples
above, the selection event was some part of the normal bundle of such a polyhedral convex set. Not all problems can be described by such polyhedral sets. In
this section we give some examples with non-trivial curvature. As above, many
of the interesting selection events are still parts of the normal bundle.
5.1. Projection onto a convex set
The problem of projection onto a closed convex set K is canonical example of
a convex program. Hence, it seems to be a natural place to start to introduce
curvature into our discussion. (Indeed two variants of the LASSO above can be
expressed in terms of projection onto a convex set.) Let
1
kβ − sk22
2
(
0 β∈K
P(β) =
∞ otherwise.

ℓ(β; s) =

In this case, with ǫ = 0 the randomized program takes the form
1
minimize kβ − sk22 − ω T β.
β∈K
2
with solution
β̂(s, ω) = PK (s + ω),

where PK is the Euclidean metric projection onto the set K (i.e. the projection
with metric induced by the ℓ2 norm).
The subgradients to our penalty are
∂P(β) = Nβ (K),
the normal cone of K at β.
In this case, the subgradient equations take the form
s + ω = β + z,

β ∈ K, z ∈ Nβ (K).

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

36

The normal cone Nβ (K) can be parameterized by its unit vectors
S(Nβ (K)) = {u ∈ Nβ (K) : kuk2 = 1}
and their length r ∈ [0, ∞).
We can therefore draw from the joint distribution of (s, ω) with a density
proportional to
f (s) · g(β + r · u − s) · |Jψ (β, r, u)|,
where4

ψ(s, β, r, z) = (s, β + r · u − s, β, β − s, r · u)

β ∈ K, u ∈ S(Nβ (K)), r ≥ 0.

Integrating out s yields the density
(f ∗ g)(β + r · u) · |Jψ (β, r, u)|.
where f ∗ g is the density of the convolution F ∗ G.
The Jacobian here is recognizable as that in the Weyl-Steiner volume of
tubes formula Weyl (1939); Hotelling (1939); Johnstone and Siegmund (1989);
Takemura and Kuriki (2002); Sun (1993); Adler and Taylor (2007)
Jψ (β, r, u) = det (I + r · C−u ) = det (I + C−r·u ) ,

(65)

where Cη is the curvature matrix of ∂K at β in the direction η normal to ∂K.
Suppose now K is compact and that F ∗ G is the uniform distribution on
K ⊕ δB2 . This can be achieved by taking F = δ0 , say and G the uniform
distribution itself. In this case, supp(G) = K ⊕ δB2 and we arrive at WeylSteiner’s tube formula
Z
1
Jψ (β, r, u) dβ dr du
1=
Hp (K ⊕ δB2 ) (β,r,u):β+r·u∈K⊕δB2

Z δ Z
=
Jψ (β, r, u) dβ du dr.
0

The formula above must be interpreted as a sum over pieces or strata of K
of different dimensions so that dβ properly refers to Hausdorff measure of differing dimensions and the matrix in (65) is also of differing dimensions. We
refer the readers to Takemura and Kuriki (2002); Adler and Taylor (2007) for
further details and points of entry into the volume-of-tubes literature. When
K is a smooth body, i.e. its boundary is a differentiable (p − 1)-dimensional
hypersurface then there is only stratum and the tube formula reads
Z δZ
Hp (K ⊕ δB2 ) =
Jψ (β, r, ηβ ) dβ dr,
0

∂K

where ηβ is the outward pointing unit normal vector field.
4 This parameterization is slightly different than our usual one, in that it constructs s + ω
instead of ω, but the difference is not important.

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

37

When F ∗ G is not the uniform distribution on K ⊕ δB2 we might instead try
to compute
Z
(f ∗ g)(β + r · u) · Jψ (β, r, u) dβ dr du
(F ∗ G)(K ⊕ δB2 ) =
(β,r,u):β+r·u∈K⊕δB2

=

Z

δ

0

Z


(f ∗ g)(β + r · u) · Jψ (β, r, u) dβ du dr.

This problem is considered in some generality when f ∗ g is a smooth density
in Chapter 10 of Adler and Taylor (2007). When F ∗ G = N (0, Ip ), their probability, expanded in a power series in δ plays an important role in the Gaussian
Kinematic Formula Taylor (2006); Taylor and Adler (2009).
Unlike in the Weyl-Steiner formula, while the first steps are similar and are
what Weyl said was something any student of calculus could do Weyl (1939),
the above probability are not Riemannian invariants of K but depend on how
the law F ∗ G relates to K.
While the connections to curvature measures and volume-of-tubes formulae
are enlightening, we feel concrete examples are also very important. In this
section, we consider the group LASSO as a canonical example of a statistical
learning problem with curvature.
5.2. Forward stepwise with groups of variables: Kac-Rice with
groups
This example is a generalization of the first-step of forward stepwise when we
allow groups of variables to enter, and is a second example of the Kac-Rice test
described above. Non-randomized approaches to this problem are considered in
Loftus and Taylor (2014); Yang et al. (2016).
The randomized version of the Kac-Rice objective with a group LASSO
penalty determined by a partition {1, 2, . . . , p} = ∪g∈G g can be written as
maximize
η T (X T y + ω) − IK (η), y × ω ∈ F × H
p
η∈R

where
K = {η ∈ Rp :

X

g∈G

(66)

λg kηg k2 ≤ 1}

and IK is the usual characteristic function of set K. In this section we have
denoted the law of ω by H with density h to distinguish the density from the
group index g. In this case, the optimal η = η ∗ is given by

T

 Xg y + ω g
if g = g ∗
ηg∗ = λg kXgT y + ωg k2

0
otherwise,
where

g ∗ = argmax
g∈G

1
kXgT y + ωg k2 .
λg

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

38

Conditioning on the selection event


Êg∗ = (y, ω) : g ∗ = argmaxkXgT y + ωg k2 /λg ,
g∈G

the subgradient equation leads to the reconstruction map


zg ∗
T
∗
∗
∗
ω = φg (y, c, zg , (zg )g6=g ) = −X y + c
,
(zg )g6=g∗
subject to constraints
c > 0, kzg∗ k2 = λg∗ , kzg k2 ≤ λg for g 6= g ∗ .
Thus we sample (y, c, zg∗ , (zg )g6=g∗ ) from a density proportional to
f (y) · h(cz − X T y) · c|g

∗

|−1

and supported on
∗

supp(F ) × R+ × λg∗ S(R|g | ) ×

Y

λg B2 (R|g| ),

g∈G,g6=g∗
∗

∗

where λg∗ S(R|g | ) is a sphere with radius λg∗ in R|g | and λg B2 (R|g| ) is an ℓ2
ball in R|g| . The Jacobian here, is the determinant of the derivative of the map

zg∗ c 7→ −XgT∗ y + czg∗ ,

∗
∗
∗
which equals det cVg∗ zg∗ = λg∗ c|g |−1 , where Vg∗ ∈ R|g | × R|g |−1 be an
orthonormal basis completion of zg∗ .
5.3. Group LASSO
As mentioned above, the LASSO (23) can be expressed in terms of metric projection, as can many problems when ℓ is squared-error loss and P is either a
(semi-)norm in or a constraint on a seminorm. When the seminorm is polyhedral the Jacobian in our parameterization typically depends only on the design
matrix X. If the unit ball of the seminorm is not polyhedral, the Jacobian of the
parameterization involves a curvature term which depends on the design matrix
X.
We consider the group LASSO Ming and Lin (2005). The group LASSO
norm, defined by a partition
∪g∈G g = {1, . . . , p}
and weights (λg )g∈G is defined as
P(β) = PG,(λg )g∈G (β) =

X

g∈G

λg kβg k2 ,

(67)

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

where βg = β[g] are the coefficients in the g-th group.
The group LASSO problem is defined by
X
1
2
ky
−
Xβk
+
λg kβg k2 .
minimize
2
β∈Rp
2

39

(68)

g∈G

Losses other than squared-error can of course be used with this penalty – it is the
penalty that is the group LASSO. The square-root group LASSO Bunea et al.
(2013) replaces the squared-error loss above with the ℓ2 loss ky − Xβk2 .
The randomized version of the group LASSO is
X
1
ǫ
minimize
ky − Xβk22 − ω T β + kβk22 +
λg kβg k2
(69)
p
β∈R
2
2
g∈G

where ω ∼ H with density h.
Having solved (69) with solution β̂(X, y, ω), we define the set of active groups
as
n
o
Ê(X, y, ω) = g : β̂g (X, y, ω) 6= 0 .
The canonical event of interest here is
n
o
(X, y, ω) : Ê(X, y, ω) = E .

(70)

On this event, the subgradient equation gives the randomization reconstruction map
ω = φE (y, (βg )g∈E , (zh )h∈−E )




(λg βg /kβg k2 )g∈E
(βg )g∈E
,
− XT y +
= (X T X + ǫI)
(zh )h∈−E
0

(71)

with the restriction on the inactive subgradients as
kzh k2 ≤ λh , h ∈ −E.
We may choose to parameterize
{(y, ω, β, α, z) ∈ S̄0F (G, ǫ, ℓ, P) : βh = 0 ∀h ∈ −E}
in terms of (y, (βg )g∈E , (zh )h∈−E ) as
ψE (y, (βg )g∈E , (zh )h∈−E )


 



(βg )g∈E
(βg )g∈E
(λg βg /kβg k2 )g∈E
= y, (X T X + ǫI)
,
− XT y +
,
0
0
(zh )h∈−E


 

(βg )g∈E
(λg βg /kβg k2 )g∈E
− XT y − XT
,
0
(zh )h∈−E
with the restriction on the inactive subgradient
Y
(zh )h∈−E ∈
λh B2 (R|h| ),

(72)

h∈−E

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

40

where B2 (R|h| ) is the ℓ2 ball in R|h| . The density to sample from thus, is proportional to
f (y) · g(φE (y, (βg )g∈E , (zh )h∈−E )) · D((βg )g∈E ,(zh )h∈−E ) φE ,

(73)

M = ∪g∈E g, the set of active predictors, and their
restricted to (72). Denoting
P
cardinality as |M | = g∈E |g|, where |g| is size of each group, the non-trivial
Jacobian term is precisely given by
D((βg )g∈E ,(zh )h∈−E ) φE
 T
XM XM + ǫI + D1
= det
T
X−M
XM

0|M|×(p−|M|)
Ip−|M|

T
= det(XM
XM + ǫI + D1 ),



(74)

where D1 is a block diagonal matrices of the form


!!
βg βgT
λg
 ∈ R|M|×|M| .
D1 = diag 
I−
kβg k2
kβg k22
g∈E

Remark 14. For groups that are singletons (|g| = 1 for all g ∈ G), the matrix
D1 is zero as βg /kβg k2 ∈ {1, −1}.
5.4. Conditioning on active directions in group LASSO
If we choose to condition additionally on the active directions
(
)
βg
ug =
, g∈E ,
kβg k2
writing
βg = γg · ug , γg > 0,
g ∈ E;
kzh k2 ≤ λh ,
h ∈ −E,
we may express the reconstruction in (71) in terms of (y, (γg )g∈E , (zh )h∈−E ) as
φ(E,(ug )g∈E ) (y, (γg )g∈E , (zh )h∈−E )




(λg ug )g∈E
(γg ug )g∈E
T
T
−X y+
= (X X + ǫI)
(zh )h∈−E
0
restricted to
(y, (γg )g∈E , (zh )h∈−E ) ∈ Rn × (0, ∞)|E| ×

Y

λh B2 (R|h| ).

h∈−E

In this case the density we sample from is proportional to
f (y) · g(φ(E,(ug )g∈E ) (y, (γg )g∈E , (zh )h∈−E )) · D((γg )g∈E ,(zh )h∈−E ) φ(E,(ug )g∈E ) ,
(75)
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

41

restricted to support above. The Jacobian of the map in this case can be computed as polynomial in (γg )g∈E by taking a geometric approach outlined below.
Hence, conditioning further on the active directions allows us to sample from a
log-concave density.
To see an explicit computation of the Jacobian, let Vg ∈ R|g|×(|g|−1) , g ∈ E,
denote a matrix whose columns are orthonormal and orthogonal to ug and set
V = diag ((Vg )g∈E ) ∈ R|M|×(|M|−|E|) , U = diag((ug )g∈E ) ∈ R|M|×|E| .
Finally set
Γ = diag((γg I|g| )g∈E ) ∈ R|M|×|M| .
The Jacobian is the derivative of the map
(U, γ) 7→ ∇ℓ̃(ΓU ) + ΛU,
where
T
∇ℓ̃(ΓU ) = (XM
XM + ǫI|M| )ΓU,

and
Λ = diag((λg Ig )g∈E ) ∈ R|M|×|M| .
Differentiating first with respect to V (i.e. tangent to U ) and then γ, the
derivative can be written as


 
(∇2 ℓ̃(ΓU )Γ + Λ)V ∇2 ℓ(ΓU )U = ∇2 ℓ̃(ΓU )Γ + ΛV V T V Γ−1 U .
T
Since ℓ̃ is quadratic, Q = ∇2 ℓ(ΓU ) = XM
XM + ǫI|M| does not involve ΓU , the
determinant of the above expression is equivalent to

 T 
 5
V
−1
(QΓ + Λ)V QU .
det
Q
UT

Letting

Γ− = diag((γg I|g|−1 )g∈E ), Γ̄ =

 −
Γ
0

the Jacobian reduces to computing
 T −1
 
ΛV
−1 V Q
det Γ̄ I|M| + Γ̄
Z T Q−1 ΛV

= det Γ− + V T Q−1 ΛV .

0
I|E|

0
0





,

(76)

Conditioning on U is now straightforward as V will then be fixed in the above
density.
5 Since

V

U

T

we condition on the active directions, U and V are fixed matrices hence
Q−1 is fixed as well.

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

42

5.5. General geometric approach
More generally, one can take a geometric approach to obtain the Jacobian, even
when we do not condition on the active directions.
For the group LASSO, in computing the Jacobian it is convenient to assume
that ℓ is somewhat arbitrary but twice-differentiable and to introduce a new
Riemannian metric on the fiber of optimization variables ((γg )g∈E , (zg )g∈E )
over the data point s. The convenience is that the structure of the Jacobian
is similar for many problems in which the penalty is a seminorm appearing
in Lagrange form. We expect a similar formula to hold when the penalty is a
seminorm appearing in bound form. In the interest of space, we do not pursue
this here.
The Riemmanian metric we use is the pull-back of the following metric on
Rp :
−1
Zβ .
(77)
hVβ , Zβ i(s,β) = VβT ∇2 ℓ(β; s) + ǫI

With ǫ = 0 and ℓ the log-likelihood of an exponential family, in local coordinates
the above metric is inverse of the observed information metric Efron and Hinkley
(1978) evaluated at the point estimate β. It is perturbed by ǫ times the Euclidean
metric above as we have modified our convex program in randomizing it.
For the group LASSO, we then now differentiate the parameterization in
coordinates γE , (zg )g∈G . The main reason for introducing this metric is that the
image tangent vectors in the face


!
Y
Y
|g|
|h|
F =
λg · S(R ) ×
λh B2 (R )
g∈E

h∈−E

P
remain orthogonal to the normal vectors g∈E γg zg . The next theorem derives
the Jacobian when curvature component is present in the geometry.
Theorem 2 (Jacobian meta-theorem). Consider the randomized version of the
problem
minimize
ℓ(β; S) + kβk
p
β∈R

for some F -a.s. twice-differentiable convex ℓ and some seminorm
kβk = sup β T z.
z∈K

Let
• F be a smooth face of K;
• (ηi )i∈I be an orthonormal (in the Euclidean metric) frame field normal to
F;
• (Vj )j∈A be an orthonormal (in the Euclidean metric) frame field tangent
to F ,

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

43

and let ψ denote the parameterization
ψ(s, z, γ) = ∇ℓ(β(z, γ); s) + ǫ · β(z, γ) + z
P
restricted to N (F ) with β(z, γ) = i∈I γi ηi . Then,


∂
ψ
= δil
ηi ,
∂γl
(s,β)


∂
Vj ,
ψ
=0
∂γi
(s,β)

(78)

hVj , Vk ψi(s,β) = S−β (Vj , Vk ) + hVj , Vk i(s,β) ,
where Sη is the shape operator of F as it sits in Rp . Finally,
Jψ(s, z, γ) = det(∇2 ℓ(β; s) + ǫ · I) · det(G(s, z, γ) + C−β (z, γ))
= det(I + G(s, z, γ)−1 C−β (z, γ))/ det(H(s, z, γ)),

where

(79)

G(s, z, γ)jk = hVj,z , Vk,z i(s,β)

⊥
⊥
ηl,z i(s,β) ,
ηi,z , P(s,β)
H(s, z, γ)il = hP(s,β)

⊥
where P(s,β)
is projection orthogonal to the tangent space of F in the metric
(77) and C−β (z, γ) is the Euclidean curvature matrix of F in the basis Vi,z , 1 ≤
i ≤ dim(F ).

Proof. The reconstruction map reads
ω = ∇ℓ
P

X

γi ηi ; s

i

!

+ǫ

X

γi ηi

i

!

+z

γi ηi is normal to K at z ∈ F . We now compute
!
!
X
∂
2
γi ηi ; s + ǫ · I ηi
ψ(s, z, γ) = ∇ ℓ
∂γi
i
!
!
!
X
X
2
γi ηi + Vj .
γi ηi ; s + ǫ · I Vj
Vj ψ(s, z, γ) = ∇ ℓ

with the constraint

i

i

i

This is enough to establish (78). Now, the determinant we want to compute is
the determinant of the matrix


∂
ψ
V ψ ∂γ
where V ψ is the matrix whose columns are Vj ψ, and
∂
columns are ∂γ
ψ.
i

∂
∂γ ψ

is the matrix whose

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

44

P
Factoring out ∇2 ℓ ( i γi ηi ; s) + ǫI and then multiplying by the orthogonal
matrix
 T
V
γT
shows that the determinant we want to compute is det(∇2 ℓ(β; s) + ǫI) times
the determinant of the matrix


S−β (Vj , Vk ) + hVj , Vk i(s,β) 0
.
...
I
Relation (79) now follows, with the second display simply a familiar formula for
the determinant expressed in block form
det(∇2 ℓ(β; s) + ǫ · I)−1 = det(G(s, z, γ)) · det(H(s, z, γ)).
When the loss is quadratic, or we use a quadratic approximation for the loss,
we see that if we condition on z, then the term det(H(s, z, γ)) is constant and
only the first term must be computed.
Remark 15. We note that while F may not be convex, when ℓ is quadratic
(or a quadratic approximation is used for inference) and we condition on z the
distribution for inference is log-concave in γ. More generally, if f is logconcave
in s as well, then the relevant distribution of (s, γ)|z is jointly logconcave in
(s, γ).
This general theorem can be applied to compute the Jacobian of the group
LASSO with a parametrization in terms of the active and inactive directions
and the magnitudes of active coefficients. We refer the readers to Appendix D
for a rederivation of this Jacobian using the above geometric perspective.
6. Multiple views of the data
Often, an analyst might try fitting several models to a data set. These models
might have different numbers of parameters, and different objective functions.
Examples include fitting a regularization path Friedman et al. (2010) or stability
selection Meinshausen and Bühlmann (2010).
Nevertheless, if each model is an instance of (33), then there is a straightforward procedure to construct a selective sampler.
Given a collection (Gi , ǫi , ℓi , Pi )i∈I , the appropriate density is proportional
to
!
Y
f (s) ·
gi (ǫi βi + αi + zi ) · Jψi (s, βi , αi , zi )
(80)
i∈I

supported on some subset of
G

s∈supp(F )

Y

i∈I

pi

!

{(βi , αi , zi ) : βi ∈ R , αi ∈ ∂ℓi (βi ; s), zi ∈ ∂Pi (βi )} .

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

45

This set has the form of a bundle with base space supp(F ) and fibers as described
above. Typically, each ψi will be restricted to some canonical event Bi . For
example, in something similar to stability selection, each Bi might identify the
set of variables and signs chosen by the i-th randomized LASSO program.
We see that given s, the variables in the fiber are independent. Hence, these
can be sampled in parallel. For instance, on different machines initiated with
the same random seed, we can sample s IID from density f , then run Gibbs or
other samplers to sample (βi , αi , zi ).
6.1. Multiple steps of forward stepwise
As an example, we consider taking K steps of forward stepwise. The Kac-Rice
test can be extended to K steps of forward stepwise, where the selection event
is characterized by a sequence of indices with the corresponding signs that constitute the active set at step K. To make it more explicit, say we consider a sequence of active variables (j1 , . . . , jK ) with corresponding signs (s1 , s2 , . . . , sK ),
in the order in which the randomized variables enter the model. The k-th optimization in the set of K optimizations can be written as
T
⊥
maximize η T (X−A
PA
y + ωk ) − IKk (η), where y × ωk ∼ F × Gk , (81)
k−1
k−1
η∈Rp−k+1

Ak = {j1 , j2 , . . . , jk } is the active set including the k-th step, X−Ak are the
columns of X except for the ones corresponding to the current active set Ak ,
the characteristic function
(
0 if η ∈ Kk
IKk (η) =
∞ otherwise
and
Kk = {η ∈ Rp−k+1 : kηk1 ≤ 1}.

Here, {ωk }K
k=1 is a sequence of independent randomization variables with ωk
coming from a given distribution Gk in Rp−k+1 and corresponding density gk .
⊥
The projection PAk is onto XAk and PA
is the residual after this projection.
k
The selection event after K steps can be written as
Ê{(sk ,jk )}K
=
k=1



K
Y

n
⊥
∈
R
×
Rp−k+1 : sign(XjTk PA
y + ωk,jk ) = sk ,
y, {ωk }K
k=1
k−1
k=1

⊥
sk (XjTk PA
y
k−1

+ ωk,jk ) ≥ max
c

j∈Ak−1

⊥
|XjT PA
y
k−1


+ ωk,j |, k = 1, . . . , K .

The randomization reconstruction map for the k-th step, from the subgradient
equation is given by
T
⊥
φk (y, zk ) = −X−A
PA
y + zk ,
k−1
k−1

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

46

where, sub-differential zk ∈ Rp−k+1 from the k-th step is restricted to
zk ∈ ∂IKk (ηk∗ ),
and ηk∗ ∈ Rn−k+1 is the optimal solution for the k-th optimization, as stated in
(81). More explicitly, the normal cone is given by
∂IKk (ηk∗ ) = {c · u : u ∈ Rp−k+1 , ujk = sk , |uj | ≤ 1 ∀j ∈ Ack−1 , c > 0}.

(82)

The sampler density is thus proportional to
f (y) ·

K
Y

k=1



T
⊥
P
y
,
gk zk − X−A
Ak−1
k−1

(83)

supported on
(y, z1 , . . . , zk ) ∈ Rn ×

K
Y

∂IKk (ηk∗ ).

k=1

For logistic regression, one can replace the T statistics above with the score
statistics as described in Remark 13.
Further examples of algorithms that choose a variable based on multiple views
of the data are presented in Appendix E.
7. Selective sampling via projected Langevin
In the examples above, we sample from a joint density of data S ∈ Rn and optimization variables T ∈ Rp , from a probability space on (S, T ), achieved via a
reparametrization. The selective sampler now samples from a pull back measure
on a transformed probability space under a reconstruction map for randomization ω. This idea implemented on randomized convex programs typically allow
us to sample from a joint density that is supported on a relatively simpler region.
On a canonical selection event B(E,zE ) , defined by the active variables E and
their signs/directions zE , the joint sampling density is
h(s, t) ∝ f (s) · g(φ(E,zE ) (s, t)) · Dt φ(E,zE ) (s, t) · 1D(E,zE ) (s, t),

(84)

supported on constraint set D(E,zE ) ⊂ Rn × Rp , under reconstruction map
φ(E,zE ) .
When we consider
f (s) ∝ exp(−f˜(s)) and g(ω) ∝ exp(−g̃(ω)),
the negative of logarithm of the sampling density is proportional to
h̃(s, φ(E,zE ) (s, t)) = f˜(s)) + g̃(φ(E,zE ) (s, t)) − log Dt φ(E,SE ) (s, t)

(85)

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

47

supported on set D(E,zE ) .
We sample from the target density h̃ using updates from a projected Langevin
random walk.
Projected Langevin Updates
Iterative update [k+1]: The (k + 1) update based on previous update (s(k) , t(k) ) is given by
!
!!
 (k+1) 
 (k) 
(k)
p
∇s h̃(s(k) , φ(E,zE ) (s(k) , t(k) ))
s
s
ξ1
+ 2η
=P
−η
(k)
t(k+1)
t(k)
∇t h̃(s(k) , φ(E,zE ) (s(k) , t(k) ))
ξ2
(k)

(k)

for step-size η, ξ1 × ξ2 ∼ N (0, In ) × N (0, Ip ) and P is projection onto set D(E,zE ) .
Computing the (k + 1) update based on previous update (s(k) , t(k) ) thus involves two steps:
1: Computing gradient of the negative of log density w.r.t. (s, t), that is
!
∇s h̃(s(k) , φ(E,zE ) (s(k) , t(k) ))
∇t h̃(s(k) , φ(E,zE ) (s(k) , t(k) ))
2: Compute projection of update from a noisy version of gradient descent of the log density
!
!
 (k) 
(k)
p
∇s h̃(s(k) , φ(E,zE ) (s(k) , t(k) ))
s
ξ1
+
−
η
2η
(k)
t(k)
∇t h̃(s(k) , φ(E,zE ) (s(k) , t(k) ))
ξ2
onto constraint set D(E,zE ) .

When h̃(s, φ(E,zE ) (s, t)) is a convex function on restriction D(E,zE ) and satisfies smoothness properties
∇h̃(s1 , ω1 ) − ∇h̃(s2 , ω2 ) ≤ β |(s1 , ω1 ) − (s2 , ω2 )|
∇h̃(s1 , ω1 ) ≤ C,
for all (s1 , ω1 ), (s2 , ω2 ) ∈ D(E,zE ) , and the support D(E,zE ) is convex with nonempty interior and contained in a finite Euclidean ball, the projected Langevin
sampler indeed converges to the target density in (84) as proved in Bubeck et al.
(2015).
Remark 16. For most examples, data S is sampled from a Gaussian density,
in which case, ∇f˜ is unbounded. One should be able to remove this condition
by considering a restriction of f˜ to a bounded set of probability close to 1 under
Gaussian density f .
Remark 17. The choices for randomizations typically include a Gaussian,
Laplacian or Logistic distribution.6 Under a Logistic density, g̃ is seen to satisfy
the smoothness conditions. Qualitatively, our samples do not change much with
any of these choices of randomization.
6 Tian and Taylor (2015) use heavy-tailed distributions such as Laplace or logistic as choices
for randomization.

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

48

Remark 18. Typically, the support set D(E,zE ) is convex, but not bounded.
We could again remove the boundedness condition on D(E,zE ) , by considering a
compact subset of probability close to 1 under sampling density h.
7.1. Examples
We revisit few examples like the LASSO with fixed and random design, the KacRice with forward stepwise, the ℓ1 -penalized Logistic to implement the projected
Langevin sampler. We offer inference on the parameters in the selected model
E conditional on the canonical selection event B(E,zE ) of selecting active set
E and their signs/ directions, obtained upon solving the corresponding randomized programs. We base our tests on the randomized pivot, developed in
Tian and Taylor (2015) for selective inference in a randomized setting 7 .
Example 7.1. LASSO with fixed design (Section 4.1): The response y is modeled as N (Xb, σ 2 ) with |supp(b)| = s and σ 2 = 1. We test the null hypothesis
H0 : bE,j = 0, bE,j denoting the j-th coefficient in the selected model E, conditioning on the selection event B(E,zE ) and sufficient statistics PE\j y, corresponding to nuisance parameters as discussed in A. Thus, the sampling density
is proportional to (44) on the support
{(y, βE , u−E ) : PE\j y = Yobs , diag(zE )βE > 0, ku−E k∞ ≤ 1}.
The updates from projected Langevin sampler can be written as

√
  (k)

(k)
(k)
(k)
y − ηy (k) /σ 2 − η∇y g̃(φ(E,zE ) (y (k) , βE , u−E )) + 2ηξ1
ỹ (k+1)
√

(k)
(k)
(k)
(k)
 (k+1)  
,
P1 βE − η∇βE g̃(φ(E,zE ) (y (k) , βE , u−E )) + 2ηξ2
=
βE




√
(k+1)
(k)
(k)
(k)
(k)
(k)
u−E
P2 u−E − η∇u−E g̃(φ(E,zE ) (y , βE , u−E )) + 2ηξ3
(k)

(k)

(k)

where (ξ1 , ξ2 , ξ3 ) ∼ N (0, Ip ) × N (0, I|E| ) × N (0, Ip−|E| ) and independent
|E|
of everything else and P1 is projection onto the orthant RzE and P2 is the
p−|E|
projection onto the cube [−1, 1]
. The conditioning on nuisance statistic is
implemented by fixing PE\j y with the update being
⊥
y (k+1) = Yobs + PE\j
ỹ (k+1) ,

where PE\j is the projection onto the column space of XE\j .
Example 7.2. ℓ1 -penalized logistic & LASSO with random design matrix (Section 4.3): In ℓ1 -penalized logistic the response y is modeled as i.i.d. Bernoulli(π(xTi b)),
with |supp(b)| = s, with the same null hypothesis as above. In the case of LASSO
7 To construct the pivots from Tian and Taylor (2015) we can either integrate over the null
statistic or condition on it. In the case of LASSO or ℓ1 -penalized logistic with random design,
the null statistic is T−E defined in (50). In the case of LASSO with fixed X, the null statistic
T y. Here we report the results where we integrate over the null statistic (conditioning
is X−E
on the null is easier computationally since it reduces the size of the sampling space).

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

49

with random design, the response is modeled as y ∼ N (Xb, σ 2 In ) with σ 2 = 1.
If we denote as β̄E the unpenalized MLE using only selected covariates XE in
the respective problems, then we can describe the inference for both problems at
once. We sample from the selective density in (53), supported on
{(t, βE , u−E ) : PE\j t = Tobs , diag(zE )βE > 0, ku−E k∞ ≤ 1},
conditioning again on selection event B(E,zE ) and nuisance statistics. Here,


Σ̃E,E 0
PE\j is the projection onto the column space of
and the nuisance
0
0
statistic equals
PE\j T = Σ̃E,E β̄E
where ΣE,E is the covariance of β̄E (we use a bootstrap estimate of this covariΣE,E eE,j eT
E,j
E,j ΣE,E eE,j

ance in the shown results) and Σ̃E,E = ΣE,E − eT

. The updates in the

projected Langevin sampler are

√
  (k)

(k)
(k)
(k)
t −ηΣ−1 t(k) − η∇t g̃(φ(E,zE ) (t(k) , βE , u−E )) + 2ηξ1
t̃(k+1)
√

(k)
(k)
(k)
(k)
 (k+1)  
,
P β − η∇βE g̃(φ(E,zE ) (t(k) , βE , u−E )) + 2ηξ2
=
βE
 1 E


√
(k+1)
(k)
(k)
(k)
(k)
u−E
P2 u−E − η∇u−E g̃(φ(E,zE ) (t(k) , βE , u−E )) + 2ηξ3
(k)

(k)

(k)

where (ξ1 , ξ2 , ξ3 ) ∼ N (0, Ip ) × N (0, I|E| ) × N (0, Ip−|E| ) and independent
|E|
of everything else and P1 is projection onto the orthant RzE and P2 is the
p−|E|
projection onto the cube [−1, 1]
. The conditioning as usual is implemented
by keeping PE\j t fixed and updating
⊥
t̃(k+1) .
t(k+1) = Tobs + PE\j

Example 7.3. Forward stepwise for Kac-Rice (Section 6.1): We perform K
steps of randomized forward stepwise as in (81) with the response generated as
y ∼ N (Xb, σ 2 In ) and |supp(b)| = s. The null hypothesis is H0 : bjK = 0, where
jk is the predictor chosen by the algorithm in the kth step, k = 1, . . . , K. We
sample from the density (53) with s = K − 1 supported on

(y, z1 , . . . , zK ) : PAK−1 y = Yobs , zk ∈ ∂IKk (ηk∗ ), k = 1, . . . , K ,

where we condition on the sufficient statistic for the nuisance parameters in
the selected model consisting of predictors XAK . The updates from projected
Langevin sampler when one computes the Kac-Rice objective conditional on the
selection event determined by forward stepwise with K steps are
P
 √


(i)
(i)
K
 (i+1)
g̃k (φk (y (i) , zk )) + 2ηξ1
y (i) − ηy (i) /σ 2 − η∇y
k=1
ỹ




√
(i)
(i)


(i) (i)
 z1 
2ηξ
P
g̃
(φ
(y
,
z
))
+
z
−
η∇


1
z1 1 1
2
1
1
 
,

=
 .. 

..
 . 


.




√
zK
(i)
(i)
(i) (i)
PK zK − η∇zK g̃K (φK (y , zK )) + 2ηξK+1
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

50

where Pk is the projection onto ∂IKk (ηk∗ ), given in (82), for all k = 1, . . . , K.8
The conditioning step is then done as above by updating y while keeping PAK−1 y
fixed.
Example 7.4. Group LASSO with fixed X (Section 5.4): The response y is
modeled as N (Xb, σ 2 ) with |supp(b)| = s and σ 2 = 1. We test the null hypothesis
H0 : bE,g = 0, bE,g denoting the coefficient corresponding to group g, g ∈ E,
conditioning on the selection event B(E,(ug )g∈E ) .9 Additionally, we condition
on the sufficient statistics for the nuisance parameters, PE\j y. The sampling
density is proportional to (75) on the support
)
(
Y
n
|E|
|h|
(y, (γg )g∈E , (zh )h∈−E ) ∈ R × (0, ∞) ×
λh B2 (R ) : PE\g y = Yobs .
h∈−E

The updates for



ỹ (k+1)
(k+1) 

 (γg )g∈E 
(k+1)
(zh )h∈−E

from projected Langevin sampler can be written as


√
(k)
(k)
(k)
2
(k)
y
−
ηy
/σ
−
η∇
g̃(φ
)
+
2ηξ
y
1

√
 
(k)
(k) 
P1 (γg )(k)
− η∇(γg )g∈E g̃(φ(k) ) + η(Tr(Dg ))g∈E + 2ηξ2  ,
g∈E




√
(k)
(k)
P2 (zh )h∈−E − η∇(zh )h∈−E g̃(φ(k) ) + 2ηξ3
where

(k)

(k)

φ(k) = φ(E,(ug )g∈E ) (y (k) , (γg )g∈E , (zh )h∈−E ),
(k)

(k)

(k)

(ξ1 , ξ2 , ξ3 ) ∼ N (0, In ) × N (0, I|E| ) × N (0, Ip−|M| ),

independent of everything else, P1 is projectionQonto the orthant (0, ∞)|E| and
P2 is the projection onto the product of balls h∈−E λh B2 (R|h| ). Here, term
(k)

(Tr(Gg ))g∈E comes from differentiating the logarithm of the Jacobian (76)
with respect to (λg )g∈E at iteration k:

−1
Dg(k) = (Γ− )(k) + V T Q−1 ΛV
Gg ,
where

Gg =



∂V T Q−1 ΛV
and (Γ− )(k) = diag (γg(k) I|g|−1 )g∈E .
∂λg

Note that Gg is a fixed matrix, thus at every iteration we just need to invert

(k)
(Γ− )(k) + V T Q−1 ΛV in order to compute Dg . The conditioning on PE\g y
is done as in the above examples.
8 This

projection is easily done using Remark 10.
includes conditioning on the active directions as well to get a log-concave density
(Section 5.4)
9 This

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

51

Acknowledgements
Jonathan Taylor was supported in part by National Science Foundation grant
DMS-1208857 and Air Force Office of Sponsored Research grant 113039. Jelena
Markovic was supported by Stanford Graduate Fellowship. Jonathan Taylor
would like to thank the Berkeley Institute of Data Science, where part of this
manuscript was written while on sabbatical in Fall 2015.
References
Adler, R. J. and Taylor, J. (2007), Random Fields and Geometry, 1 edn,
Springer.
Barber, R. F. and Candes, E. (2014), ‘Controlling the false discovery rate via
knockoffs’, arXiv:1404.5609 [math, stat] . arXiv: 1404.5609.
URL: http://arxiv.org/abs/1404.5609
Barndorff-Nielsen, O. (1983), ‘On a formula for the distribution of the maximum
likelihood estimator’, Biometrika 70(2), 343–365.
URL: http://biomet.oxfordjournals.org/content/70/2/343
Becker, S., Candes, E. and Grant, M. (2010), Templates for convex cone problems with applications to sparse signal recovery, Technical report, Stanford
University.
URL: http://arxiv.org/abs/1009.2065
Belloni, A., Chernozhukov, V. and Wang, L. (2014), ‘Pivotal estimation via
square-root Lasso in nonparametric regression’, The Annals of Statistics
42(2), 757–788. arXiv: 1105.1475.
URL: http://arxiv.org/abs/1105.1475
Benjamini, Y. (2010), ‘Simultaneous and selective inference: Current successes and future challenges’, Biometrical Journal. Biometrische Zeitschrift
52(6), 708–721.
Berk, R., Brown, L., Buja, A., Zhang, K. and Zhao, L. (2013), ‘Valid postselection inference’, The Annals of Statistics 41(2), 802–837.
URL: http://projecteuclid.org/euclid.aos/1369836961
Besag, J. and Clifford, P. (1989), ‘Generalized monte carlo significance tests’,
Biometrika 76(4), 633–642.
Bien, J., Taylor, J. and Tibshirani, R. (2013), ‘A lasso for hierarchical interactions’, Annals of statistics 41(3), 1111.
Bubeck, S., Eldan, R. and Lehec, J. (2015), ‘Sampling from a logconcave distribution with Projected Langevin Monte Carlo’, arXiv preprint
arXiv:1507.02564 .
URL: http://arxiv.org/abs/1507.02564
Bunea, F., Lederer, J. and She, Y. (2013), ‘The Group Square-Root Lasso:
Theoretical Properties and Fast Algorithms’, arXiv preprint arXiv:1302.0261
. arXiv: 1302.0261.
URL: http://arxiv.org/abs/1302.0261
Chan, H. P. and Walther, G. (2011), ‘Detection with the scan and the average
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

52

likelihood ratio’, arXiv preprint arXiv:1107.4344 .
URL: https://arxiv.org/abs/1107.4344
Chen, S., Donoho, D. and Saunders, M. (1998), ‘Atomic decomposition for basis
pursuit’, SIAM Journal on Scientific Computing 20(1), 33–61.
URL: http://epubs.siam.org/doi/abs/10.1137/S1064827596304010
Choi, Y., Taylor, J. and Tibshirani, R. (2014), ‘Selecting the number of principal
components: estimation of the true rank of a noisy matrix’, arXiv preprint
arXiv:1410.8260 . arXiv: 1410.8260.
URL: http://arxiv.org/abs/1410.8260
Cox, D. (1975), ‘A note on data-splitting for the evaluation of significance levels’,
Biometrika 62(2), 441–444.
URL: http://biomet.oxfordjournals.org/content/62/2/441.full.pdf
Diaconis, P., Holmes, S. and Shahshahani, M. (2013), Sampling from a Manifold,
in ‘Institute of Mathematical Statistics Collections’, Institute of Mathematical
Statistics, Beachwood, Ohio, USA, pp. 102–125.
URL: http://projecteuclid.org/euclid.imsc/1379942050
Diaconis, P. and Sturmfels, B. (1998), ‘Algebraic algorithms for sampling from
conditional distributions’, The Annals of Statistics 26(1), 363–397.
URL: http://projecteuclid.org/euclid.aos/1030563990
Dwork, C., Feldman, V., Hardt, M., Pitassi, T., Reingold, O. and Roth, A.
(2015), ‘The reusable holdout: Preserving validity in adaptive data analysis’,
Science 349(6248), 636–638.
URL: http://www.sciencemag.org/content/349/6248/636
Efron, B. and Hinkley, D. V. (1978), ‘Assessing the Accuracy of the Maximum Likelihood Estimator: Observed Versus Expected Fisher Information’,
Biometrika 65(3), 457–482.
URL: http://www.jstor.org/stable/2335893
Efron, B. and Tibshirani, R. (1996), ‘Using specially designed exponential families for density estimation’, The Annals of Statistics 24(6), 2431–2461.
URL: http://projecteuclid.org/euclid.aos/1032181161
Federer, H. (1959), ‘Curvature measures’, Transactions of the American
Mathematical Society 93(3), 418–491.
URL:
http://www.ams.org/journals/tran/1959-093-03/S0002-9947-19590110078-1/S0002-9947-1959-0110078-1.pdf
Fithian, W., Sun, D. and Taylor, J. (2014), ‘Optimal Inference After Model
Selection’, arXiv preprint arXiv:1410.2597 . arXiv: 1410.2597.
URL: http://arxiv.org/abs/1410.2597
Fithian, W., Taylor, J., Tibshirani, R. and Tibshirani, R. (2015), ‘Selective
sequential model selection’, arXiv preprint arXiv:1512.02565 .
URL: https://arxiv.org/abs/1512.02565
Frick, K., Munk, A. and Sieling, H. (2014), ‘Multiscale change point inference’,
Journal of the Royal Statistical Society: Series B (Statistical Methodology)
76(3), 495–580.
URL: http://onlinelibrary.wiley.com/doi/10.1111/rssb.12047/suppinfo
Friedman, J., Hastie, T. and Tibshirani, R. (2008), ‘Sparse inverse covariance
estimation with the graphical lasso’, Biostatistics 9(3), 432–441.
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

53

URL: http://biostatistics.oxfordjournals.org/content/9/3/432.short
Friedman, J., Hastie, T. and Tibshirani, R. (2010), ‘Regularization paths for
generalized linear models via coordinate descent’, Journal of Statistical Software 33(1), 1–22.
URL: https://core.ac.uk/download/files/153/6287975.pdf
Hillier, G. and Armstrong, M. (1999), ‘The Density of the Maximum Likelihood
Estimator’, Econometrica 67(6), 1459–1470.
URL: http://www.jstor.org/stable/2999568
Hotelling, H. (1939), ‘Tubes and spheres in n-spaces and a class of statistical
problems’, American Journal of Mathematics 61, 440–460.
Hyun, S., G’Sell, M. and Tibshirani, R. J. (2016), ‘Exact Post-Selection Inference
for Changepoint Detection and Other Generalized Lasso Problems’, arXiv
preprint arXiv:1606.03552 .
URL: https://arxiv.org/abs/1606.03552
Johnstone, I. and Siegmund, D. (1989), ‘On Hotelling’s formula for the volume
of tubes and Naiman’s inequality’, The Annals of Statistics 17(1), 184–194.
URL: https://projecteuclid.org/euclid.aos/1176347010
Lee, J. D., Sun, D. L., Sun, Y. and Taylor, J. E. (2016), ‘Exact post-selection
inference with the lasso’, The Annals of Statistics 44(3), 907–927.
URL: http://projecteuclid.org/euclid.aos/1460381681
Lee, J. D. and Taylor, J. E. (2014), Exact Post Model Selection Inference
for Marginal Screening, in Z. Ghahramani, M. Welling, C. Cortes, N. D.
Lawrence and K. Q. Weinberger, eds, ‘Advances in Neural Information
Processing Systems 27’, Curran Associates, Inc., pp. 136–144.
URL:
http://papers.nips.cc/paper/5292-exact-post-model-selectioninference-for-marginal-screening.pdf
Loftus, J. R. and Taylor, J. E. (2014), ‘A significance test for forward stepwise
model selection’, arXiv preprint arXiv:1405.3920 .
URL: https://arxiv.org/abs/1405.3920
Meinshausen, N. and Bühlmann, P. (2006), ‘High-dimensional graphs and variable selection with the lasso’, The Annals of Statistics 34(3), 1436–1462.
URL: https://projecteuclid.org/euclid.aos/1152540754
Meinshausen, N. and Bühlmann, P. (2010), ‘Stability selection’, Journal of the
Royal Statistical Society: Series B (Statistical Methodology) 72(4), 417–473.
URL: http://onlinelibrary.wiley.com/doi/10.1111/j.
1467-9868.2010.00740.x/abstract
Meinshausen, N., Meier, L. and Bühlmann, P. (2012), ‘P-values for highdimensional regression’, Journal of the American Statistical Association .
Ming, Y. and Lin, Y. (2005), ‘Model selection and estimation in regression
with grouped variables’, Journal of the Royal Statistical Society: Series B
68(1), 49–67.
URL: http://onlinelibrary.wiley.com/doi/10.1111/j.
1467-9868.2005.00532.x/abstract
Negahban, S., Ravikumar, P., Wainwright, M. J. and Yu, B. (2010), ‘A unified
framework for high-dimensional analysis of $m$-estimators with decomposable regularizers’, arXiv:1010.2731 .
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

54

URL: http://arxiv.org/abs/1010.2731
Panigrahi, S., Taylor, J. and Weinstein, A. (2016), ‘Bayesian post-selection inference in the linear model’, arXiv preprint arXiv:1605.08824 .
URL: https://arxiv.org/abs/1605.08824
Rivera, C. and Walther, G. (2013), ‘Optimal detection of a jump in the intensity
of a poisson process or in a density with likelihood ratio statistics’, Scandinavian Journal of Statistics 40(4), 752–769.
URL: http://onlinelibrary.wiley.com/doi/10.1111/sjos.12027/full
Sabourin, J., Valdar, W. and Nobel, A. (2014), ‘A Permutation Approach for
Selecting the Penalty Parameter in Penalized Model Selection’, arXiv preprint
arXiv:1404.2007 . arXiv: 1404.2007.
URL: http://arxiv.org/abs/1404.2007
Schneider, R. (1993), Convex Bodies: the Brunn-Minkowski Theory, Vol. 44
of Encyclopedia of Mathematics and its Applications, Cambridge University
Press, Cambridge.
Sun, J. (1993), ‘Tail probabilities of the maxima of Gaussian random fields’,
The Annals of Probability 21(1), 34–71.
URL: https://projecteuclid.org/euclid.aop/1176989393
Sun, T. and Zhang, C.-H. (2011), ‘Scaled Sparse Linear Regression’, arXiv
preprint arXiv:1104.4595 . arXiv: 1104.4595.
URL: http://arxiv.org/abs/1104.4595
Takemura, A. and Kuriki, S. (2002), ‘Maximum of Gaussian field on piecewise smooth domain: Equivalence of tube method and Euler characteristic
method.’, Ann. of Appl. Prob. 12(2), 768–796.
URL: https://projecteuclid.org/euclid.aoap/1026915624
Taylor, J. E. (2006), ‘A Gaussian kinematic formula’, The Annals of Probability
34(1), 122–158.
URL: https://projecteuclid.org/euclid.aop/1140191534
Taylor, J. E. and Adler, R. J. (2009), ‘Gaussian processes, kinematic formulae
and Poincares limit’, The Annals of Probability 37(4), 1459–1482.
URL: https://projecteuclid.org/euclid.aop/1248182144
Taylor, J., Loftus, J. and Tibshirani, R. (2013), ‘Tests in adaptive regression via
the Kac-Rice formula’, The Annals of Statistics 44(2), 743–770.
URL: http://projecteuclid.org/euclid.aos/1458245734
Taylor, J. and Tibshirani, R. (2016), ‘Post-selection inference for l1-penalized
likelihood models’, arXiv preprint arXiv:1602.07358 .
URL: http://arxiv.org/abs/1602.07358
Tian, X., Bi, N. and Taylor, J. (2016), ‘MAGIC: a general, powerful and
tractable method for selective inference’, ArXiv e-prints .
URL: https://arxiv.org/abs/1607.02630
Tian, X., Loftus, J. R. and Taylor, J. E. (2015), ‘Selective inference with unknown variance via the square-root LASSO’, arXiv preprint arXiv:1504.08031
. arXiv: 1504.08031.
URL: http://arxiv.org/abs/1504.08031
Tian, X. and Taylor, J. E. (2015), ‘Selective inference with a randomized response’, arXiv preprint arXiv:1507.06739 . arXiv: 1507.06739.
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

55

URL: http://arxiv.org/abs/1507.06739
Tibshirani, R. (1996a), ‘Regression shrinkage and selection via the lasso’, Journal of the Royal Statistical Society: Series B 58(1), 267–288.
Tibshirani, R. (1996b), ‘Regression shrinkage and selection via the lasso’, Journal of the Royal Statistical Society: Series B 58(1), 267–288.
URL: http://onlinelibrary.wiley.com/doi/10.1111/j.
1467-9868.2011.00771.x/abstract
Tibshirani, R. J. (2013), ‘The Lasso Problem and Uniqueness’, Electronic Journal of Statistics 7, 1456–1490.
URL: http://projecteuclid.org/euclid.ejs/1369148600
Tibshirani, R. J. (2015), ‘A General Framework for Fast Stagewise Algorithms’,
Journal of Machine Learning Research 16, 2543–2588.
URL: http://jmlr.org/papers/v16/tibshirani15a.html
Tibshirani, R. J. and Taylor, J. (2011), ‘The solution path of the generalized
lasso’, The Annals of Statistics 39(3), 1335–1371.
URL: http://projecteuclid.org/euclid.aos/1304514656
Tibshirani, R. J. and Taylor, J. (2012), ‘Degrees of freedom in lasso problems’,
The Annals of Statistics pp. 1198–1232.
Tibshirani, R., Taylor, J., Lockhart, R. and Tibshirani, R. (2014), ‘Post-selection
adaptive inference for Least Angle Regression and the Lasso’, arXiv preprint
arXiv:1401.3889 .
URL: http://arxiv.org/abs/1401.3889
Van de Geer, S. A. (2008), ‘High-dimensional generalized linear models and the
lasso’, The Annals of Statistics pp. 614–645.
Wainwright, M. J. (2009), ‘Sharp thresholds for high-dimensional and noisy
sparsity recovery using-constrained quadratic programming (lasso)’, Information Theory, IEEE Transactions on 55(5), 2183–2202.
Weyl, H. (1939), ‘On the volume of tubes’, American Journal of Mathematics
61, 461–472.
Yang, F., Foygel Barber, R., Jain, P. and Lafferty, J. (2016), ‘Selective Inference
for Group-Sparse Linear Models’, ArXiv e-prints .
URL: https://arxiv.org/abs/1607.08211
Zhao, P. and Yu, B. (2006), ‘On model selection consistency of lasso’, The Journal of Machine Learning Research 7, 2541–2563.
Zhou, Q. (2014), ‘Monte Carlo Simulation for Lasso-Type Problems by Estimator Augmentation’, Journal of the American Statistical Association
109(508), 1495–1516. arXiv: 1401.4425.
URL: http://arxiv.org/abs/1401.4425
Zou, H. and Hastie, T. (2005), ‘Regularization and variable selection via the
elastic net’, Journal of the Royal Statistical Society: Series B 67(2), 301–320.
URL: http://onlinelibrary.wiley.com/doi/10.1111/j.
1467-9868.2005.00503.x/abstract

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

56

Appendix A: Nuisance parameters & conditional distributions
In forming selective hypothesis tests or intervals when the model M is an exponential family, it was noted in Fithian et al. (2014) that the classical approach
of conditioning on appropriate sufficient statistics can be used to eliminate nuisance parameters. In this section we describe selective samplers that can be used
to sample from such conditional distributions.
Specifically, suppose we use the randomized Lasso (41) with X considered
fixed and Ē is a subset of features that we will use to form the selected model
n
o
2
2
MĒ = N (XĒ βĒ , σĒ
) : βĒ ∈ R|Ē| , σĒ
>0 .
(86)
Often, for a given (E, sE ) observed after fitting the randomized Lasso, we will
2
choose Ē = E but this is not strictly necessary. As σĒ
is part of the indexing
2
set for MĒ , we are assuming σĒ is unknown.
Suppose we want a selective test of
H0,j|Ē : βj|Ē = 0
in the model MĒ . Standard exponential family calculations detailed in Fithian et al.
(2014) tell us that we can construct such a test in the presence of the nuisance
2
parameters (βĒ\j|Ē , σĒ
) by conditioning on the appropriate sigma-algebra:



T
σ XĒ\j
y, kyk22 = σ PĒ\j y, k(I − PĒ\j )yk22 .

For any fixed values of the sufficient statistics, say (Wobs , SSEobs ) (with
PĒ\j Wobs = Wobs ) the conditional distribution for any F ∈ MĒ is supported
on the set

Dobs = y : y = Wobs + r, PĒ\j r = 0, krk2 = SSEobs .
Under H0,j|Ē it is uniformly distributed over the above set.
We must therefore sample from the set

{(y, βE , u−E ) : y ∈ Dobs , diag(sE )βE > 0, kuk−E ≤ 1}
with a density proportional to
  


βE
sE
T
f (y) · g ǫ
.
− X (y − XE βE ) + λ
0
u−E
Appendix B: LASSO and variants
In this section, we describe two common variants of the LASSO.

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

57

B.1. LASSO in bound form
The LASSO program in bound form is defined as
minimize

β∈Rp :kβk1 ≤δ

1
ky − Xβk22
2

(87)

with its corresponding randomized version
minimize

β∈Rp :kβk1 ≤δ

1
ǫ
ky − Xβk22 + kβk22 − ω T β.
2
2

(88)

In the notation established so far
S=y
1
ℓ(β; y) = ky − Xβk22
2
(
0 kβk1 ≤ δ
P(β) =
∞ otherwise.
Typically, we will be interested in doing inference when the constraint above
is tight. In this case, for the canonical event B(E,zE ) the set S̄0F (G, ǫ, ℓ, P) can
be parametrized by

(y, βE , u−E , c) : y ∈ Rn ,
βE ∈ R|E| , kβE k1 = δ, diag(zE )βE > 0,

u−E ∈ Rp−|E| , ku−E k∞ ≤ 1,

c ∈ R, c > 0 .

Another possible reparametrization, can be stated as

(y, βE\1 , u−E , c) : y ∈ Rn ,

βE\1 ∈ R|E|−1 , sE,1 βE,1 = δ − sTE\1 βE\1 , diag(sE\1 )βE\1 > 0,

u−E ∈ Rp−|E| , ku−E k∞ ≤ 1,

c ∈ R, c > 0 ,

where βE\1 = (βE,2 , . . . , βE,|E| ) and similarly sE\1 . The above parameterization
can be expressed as

 


βE
zE
T
ψ(E,zE ) (y, βE\1 , u−E , c) = y, ǫ
,
− X (y − XE βE ) + c
0
u−E
 


βE
zE
, −X T (X − XE βE ), c
,
0
u−E
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

58

where βE is expressed in terms of βE\1 as above and the corresponding randomization reconstruction map is


 
zE
β
,
φ(E,zE ) (y, βE\1 , u−E , c) = ǫ E − X T (y − XE βE ) + c
u−E
0
where



(δ − sTE\1 βE\1 )sE,1
.
βE = 
βE\1
0

The sampling density is proportional to


f (y) · g φ(E,zE ) (y, βE\1 , u−E , c) · det(D(βE\1 ,u−E ,c) φ(E,zE ) ) ,

where the last determinant equals
−c

p−|E|

−ǫsTE\1 sE,1 + X1T (−sE,1 X1 sTE\1 + XE\1 ) sE,1
det
T
ǫI + XE\1
(−sE,1 X1 sTE\1 + XE\1 )
sE\1

!

.

B.2. Basis pursuit
The basis pursuit problem Chen et al. (1998) is defined as
minimize

β∈Rp :ky−Xβk2 ≤δ

kβk1 .

(89)

In the notation established so far
S=y
(

ℓ(β; y) =

0 ky − Xβk2 ≤ δ
∞ otherwise

P(β) = kβk1 .

The set on which the problem (89) has a solution is
n
o
y : θ̂(y) 6= 0 = {y : k(I − PC )yk2 ≤ δ} ,

where PC denotes projection onto col(X).
Its randomized version is
minimize

β∈Rp :ky−Xβk2 ≤δ

ǫ
kβk22 − ω T β + kβk1 .
2

(90)

This problem, with non-random choices of ω is considered in Becker et al. (2010).
We see that
n
o

(y, ω) : θ̄(y, ω) 6= ∅ = (y, ω) : θ̂(y) 6= ∅ .
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

59

Generically, when a solution exists, the constraint will be tight, hence our
events of interest will typically condition on ky − Xβk2 = δ. We see then

S̄0F (G, ǫ, ℓ, P) = (y, ω, β, α, z) : y ∈ supp(F ), ky − Xβk2 = δ, ω ∈ supp(G),
ω = ǫβ + α + z, α = cX T (Xβ − y), c > 0, z ∈ ∂P(β) .

by

For the canonical event B(E,zE ) , the set S̄BF (G, ǫ, ℓ, P) can be parameterized

(y, βE , u−E , c) : k(I − PC )yk2 ≤ δ,

βE ∈ R|E| , diag(zE )βE > 0, ky − XE βE k2 = δ,

u−E ∈ Rp−|E| , ku−E k∞ ≤ 1,

c ∈ R, c > 0 .

(91)

This set has the form of a bundle over a subset of Rn , with fibers that are the
product a |E| − 1 dimensional ellipse, a p − |E| dimensional cube and a half-line.
The parameterization is

 


β
zE
ψ(E,zE ) (y, βE , u−E , c) = y, ǫ E − cX T (y − XE βE ) +
,
0
u−E
 


βE
zE
.
, −cX T (y − XE βE ),
0
u−E


 
zE
βE
T
, the
− cX (y − XE βE ) +
Denoting φ(E,zE ) (βE , u−E , c) = ǫ
u−E
0
sampling density is proportional to

f (y) · g φ(E,zE ) (βE , u−E , c) · det(D(βE ,u−E ,c) φ(E,zE ) )

and restricted to (91), where det(D(βE ,u−E ,c) φ(E,zE ) ) is the Jacobian of the
parameterization in the coordinates described above, i.e. coordinates on the
product of an |E| − 1 dimensional ellipse, an p − |E| dimensional cube and a
half-line.
B.3. Square-root LASSO

Using the LASSO with an unknown noise level is somewhat of a chicken-and-egg
problem, as knowing which value of λ to choose is somewhat difficult. Suppose
that instead we use the square-root LASSO. The square-root LASSO program
Belloni et al. (2014); Sun and Zhang (2011) is defined as
minimize
ky − Xβk2 + λkβk1 .
p
β∈R

(92)

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

60

Selective inference for this program was considered in (Tian et al., 2015). This
program has the advantage that it is possible to choose a reasonable value of λ
without knowing anything about the noise level. Its randomized version is
ǫ
minimize
ky − Xβk2 + kβk22 − ω T β + λkβk1 .
p
β∈R
2

(93)

We can take the law of ω to be independent of any particular noise level, though
we probably should choose ǫ = O(n−1/2 ) if thinking of something like the pairs
model.
In the notation established so far
S = y,
ℓ(β; y) = ky − Xβk2 ,
P(β) = λkβk1 .

We are most interested in inference when y − X β̂(y, ω) 6= 0 on which ℓ is
differentiable and for the canonical event B(E,zE ) = {β : diag(zE )βE > 0, β−E =
0} we can parametrize {(y, ω, β, α, z) ∈ S̄0 (G, ǫ, ℓ, P) : y ∈ Dobs , β ∈ B(E,zE ) } as
ψ(E,zE ) (y, βE , u−E ) =



 
X T (y − XE βE )
zE
βE
,
+λ
y, ǫ
−
u−E
0
ky − XE βE k2
!

 
X T (y − XE βE )
zE
βE
.
,λ
,−
u−E
0
ky − XE βE k2

|E|

with the domain Dobs × RzE × [−1, 1]p−|E|. We therefore must sample from a
density proportional to
  



 T
X T (y − XE βE )
XE R(y, βE )XE
βE
zE
f (y)·g ǫ
−
+λ
+ ǫI
· det
0
u−E
ky − XE βE k2
ky − XE βE k2
and restricted to to the domain of ψ(E,zE ) , where
R(y, βE ) = In −

(y − XE βE )(y − XE βE )T
.
ky − XE βE k22

B.4. More complex penalties: generalized LASSO and other
quadratic programs
Clearly, not every statistical learning problem of interest can be expressed in
terms of the LASSO or forward stepwise. For example, the generalized LASSO
Tibshirani and Taylor (2011) considers a penalty of the form
P(β) = λkDβk1 .
Selective inference for the entire solution path of the generalized LASSO, analgous to the solution path of Least Angle Regression and LASSO in Tibshirani et al.
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

61

(2014) has been considered in Hyun et al. (2016). We consider a fixed value of
λ combined with randomization. For some D, the programs can be re-expressed
as LASSO problems but not all D (c.f. Tibshirani and Taylor (2011)). More
generally, we might be confronted with solving a problem of the form
minimize ℓ(β; S) + P̄(Dβ − γ).
β

(94)

That is, we consider our usual problem with
P(β) = P̄(Dβ − γ).
For example, a linearly constrained quadratic program might take the form
1
minimize β T Qβ − S T β
β
2

(95)

subject to Dβ ≤ γ with S ∼ F . In this case, P̄ is the non-positive cone constraint.
Sometimes, the Fenchel conjugate P ∗ is simple in the sense that its naturally
associated selection events are easy to parametrize. In this case, we can proceed
as before and simply consider randomized programs of the form
ǫ
1
minimize β T Qβ − (S + ω)T β + P̄(Dβ − γ) + kβk22 .
β
2
2

(96)

However, when P ∗ is complex, then our usual approach may be computationally expensive. In this section, we describe an alternative randomization scheme
that may yield simpler sampling algorithms as described in Section 7.
For concreteness, we consider the quadratic program (96). Similar calculations hold for the generalized LASSO by simply replacing one convex conjugate
with another. The issue of parameterization of selection events arises quickly
upon inspection of the subgradient equation for the randomized program (96).
It will often be of interest to condition on the set of tight constraints
o
n
Ê(S, ω) = j : DjT β̂(S, ω) = γj .
On this event, the subgradient equation reads

T
ω = Qβ − S + DE
zE
T
where zE ≥ 0 and DE β = γE . We see the sub-gradient DE
zE lies in the cone
generated by the rows of DE . If we were to parametrize this selection event, we
might write it as
ω = Qβ − S + z̄E

where z̄E ∈ cone(DE ). For arbitrary D and E, this cone may be rather complex.
Each step of the projected Langevin implementation described in Section 7
requires projection onto this set. Projection onto this set can be achieved via a
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

62

non-negative least squares problem, which would in principal need to be re-run
for every step of the sampler.
We shall instead solve the corresponding randomized dual program. Introducing variable v = Dβ − γ, the dual of (95) which solves


ǫ
1
maximize minimize β T Qβ − S T β + kβk22 + uT (Dβ − γ − v) ,
u
β,v:v≤0 2
2
subject to u ≥ 0 is equivalent to solving
1
minimize (DT u − S)T (Q + ǫI)−1 (DT u − S) + uT γ
u
2
subject to u ≥ 0 (for the generalized LASSO, the non-negative cone constraint is replaced with the constraint kuk∞ ≤ λ). The set of tight constraints
can be thought of as corresponding to non-zero u’s in this dual problem E.
When the dual problem is strongly convex this is a one-to-one correspondence,
though if it is not strongly convex there may be more than one dual solution
Tibshirani and Taylor (2012).
The randomized dual problem solves the program
1
minimize (DT u − S)T (Q + ǫI)−1 (DT u − S) + uT (γ − ω)
u
2

(97)

subject to u ≥ 0, which is dual to
1
ǫ
minimize β T Qβ − S T β + kβk22 + IK (β)
β
2
2
with
K = {β : Dβ − γ + ω ≤ 0}
and IK denoting its characteristic function. However, for the above randomized
QP that we propose to solve, the law of ω has to supported on the set


ω : ∩j β : DjT β ≤ γj − ωj 6= 0 .
to ensure feasibility of the primal problem, though the dual problem always has
a solution.

Remark 19. Generally, this is not a problem in examples like the generalized
LASSO, where P is finite everywhere. Further, in the applications of selective
inference we have in mind, a data analyst is typically solving statistical learning
programs to suggest parameters of interest. For example, Heirnet Bien et al.
(2013) is a statistical learning method that seeks hierarchical interactions that
can be expressed as a linearly constrained quadratic program. In this case, a
randomization which makes the primal problem infeasible is not catastrophic.
Denoting
Qǫ = (Q + ǫI)−1 ,
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

63

one solution to the above problem that ensures feasibility of the primal problem
is to instead consider
ǫ′
1
minimize (DT u − S)T Qǫ (DT u − S) + uT (γ − ω) + kuk22
u
2
2

(98)

subject to u ≥ 0 for ǫ′ > 0. This always yields a feasible primal
1
minimize z T (DQǫ DT + ǫ′ I)−1 z − z T (DQǫ DT + ǫ′ I)−1 (DQǫ S − γ + ω),
z≤0
2
as DQǫ DT + ǫ′ I is a positive definite matrix with added positive ǫ′ . The reconstruction map in this case is given by
 


uE
0
T
ω = φE (S, u, z−E ) = DQǫ (DE
uE − S) + ǫ′
,
(99)
+γ+
0
z−E
which allows us to sample (S, uE , z−E ) from a sampling density proportional to
f (S) · g(φE (S, uE , z−E ))
|E|

p−|E|

and supported on supp(F ) × R+ × R−
, where R+ = {x ∈ R : x ≥ 0} and
R− = {x ∈ R : x ≤ 0}.
Another possibility is to consider the dual of the randomized objective (96)
given by
1
minimize (DT u − S − ω)T Qǫ (DT u − S − ω) + uT γ,
u≥0
2
which lead to reconstruction map
T
DQǫ ω = DQǫ DE
uE − DQǫ S + γ +



0

z−E



where we have conditioned on the value of the set
Ē(S, ω) = {j : ûj (S, ω) > 0}
of being E. As noted above, the set Ē(S, ω) may not be unique even when
E(S, ω) is (c.f. Tibshirani and Taylor (2012)). In this sense, we are conditioning
on a particular value of the boundary set, the boundary set determined by the
solver we use to solve this dual problem.
Denoting the density of DQǫ ω as g̃, the sampling density of (S, uE , z−E ) is
now proportional to



0
T
f (S) · g̃ DQǫ DE uE − DQǫ S + γ +
z−E
and restricted to
z−E ≤ 0, uE ≥ 0 and z−E + γ−E ∈ col(D−E ).
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

64

⊥
To avoid enforcing the implicit constraint of PD
(z−E + γ−E ) being 0 at each
−E
step of the sampler, we can sample from an approximate density for (S, uE , z−E )
⊥
that substitutes the Dirac delta operator δ0 (PD
(z−E +γ−E )) with a smoothed
−E
version, an example being
!



2
1
0
⊥
T
·exp − ′ PD−E (z−E + γ−E )
f (S)·g̃ DQǫ DE uE − DQǫ S + γ +
z−E
2ǫ
2
|E|

p−|E|

supported on supp(F ) × R+ × R−

, a much simpler constraint set.

Appendix C: Neighborhood selection
We now consider neighborhood selection (a selective sampler for which appears
in Tian et al. (2016)). Meinshausen and Bühlmann (2006) proposed neighborhood selection with the LASSO to achieve this goal. The algorithm can be
formulated as the following optimization problem, for any node i
β̂ i,λ = argmin kxi − Xβk22 + λkβk1 ,

(100)

β∈Rp ,βi =0

where with slight abuse of notation, X ∈ Rn×p is the data matrix observed from
n i.i.d observations, and xi is the i-th column of X. Choice of λ is discussed in
Chapter 3 of Meinshausen and Bühlmann (2006). Denote
B̂ = (β̂ 1 , β̂ 2 , . . . , β̂ p ),
we propose the randomized version of (100),
B̂ =

ǫ
kX − XBk2F + λkBk1 + kBk2F − ΩB,
2
B∈Rp×p ,Bii =0 ∀i∈{1,...,p}
argmin

(101)

i.i.d

where Ω = (ω 1 , . . . , ω p ) ∼ G. This is the matrix form of (33), and the KKT
conditions are decomposable across the nodes.
Suppose for node i, E i is the active set for (100), si is the corresponding
signs, λ · ui is the subgradient corresponding to the inactive variables except
the i-th and X−i is the columns of X except the i-th column. For every node i,
the i-th coordinate of β i is held to be zero, and (100) is in fact a regression of
dimension p − 1, thus
T
i
αi = −X−i
(xi − XE i βE
(102)
i ),
and the reparametrization map,

ψ(E,zE ) (X, B − , u) = (X, α + z + ǫB − , B − , α, z),
where
α = (α1 , α2 , . . . , αp ), αi from (102),
 i
s
1 2
p
i
z = (z , z , . . . , z ), z = λ i ,
u

B − ∈ R(p−1)×p is p × p matrix without the diagonal elements.
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

65

Since ω i , 1 ≤ i ≤ p, are independent, and the Jacobian
Y
T
Jψ(E,zE ) (X, B − , u) =
det(XE
i XE i + ǫI),
i∈Γ

conditioning on (E i , si ), the distribution for selective inference has the following
density,
 i

Y  β i i 
s
T
i
g ǫ E + λ i − X−i
(xi − XE i βE
f (X) ·
i)
0
u
i∈Γ
(103)
Y
T
·
det(XE
i XE i + ǫI) .
i∈Γ

After seeing the active set E = (E 1 , . . . , E p ), it is natural to choose the
selected model (i.e. the edge set E) to be E. However, since the active set E is
not necessarily symmetric, we choose the edge set E to be E ∨ , where
E ∨ = {(i, j) | Eij = 1 or Eji = 1}.
Under this model, the distribution of X is an exponential family. More specifically, conditioning on the set of edges E ∨ ,


1
T
f (X) ∝ exp − Tr(ΘX X)
2


X
X
1
= exp −
Θij xTi xj −
Θii kxi k2  .
2
∨
(i,j)∈E

i∈Γ

Note that this is an exponential family with sufficient statistics
{xTi xj , (i, j) ∈ E ∨ , kxi k2 , i ∈ Γ}.

Therefore, the law for selective inference (103) is also an exponential family with
the same sufficient statistics. To construct the UMPU tests as in Fithian et al.
(2014) for the null hypothesis H0,ij : Θij = 0, we condition on the sufficient
statistics corresponding to the nuisance parameter and sample from the law
(103).
Appendix D: Recomputing Jacobian for group lasso
We apply the Jacobian meta theorem 2 to the group Lasso to re-derive the
Jacobian in the selective sampler density from a geometric perspective. For the
tangential frame Vi,z we can split this over groups as
(ΠTg (Vg,j,zg )1≤j≤|g|−1 )g∈E = (Vg,zg )g∈E
where Πg : Rp → Rg is projection onto the g coordinates represented by the
matrix Πg ∈ R|g|×p . The the vectors Vg,j,zg are chosen to be are orthonormal
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

66

T
z = 0. For the normal frame
within Tg (λg S(Rg )) additionally satisfying Vg,j,z
g g
we can take ηg = zg /λg and
X
β=
γg zg /λg .
g∈G

In the standard basis of Tβ Rp , the metric (77) has matrix
X T X + ǫI.
The curvature matrix C−β is 0 except on the tangent spaces coming from active
groups, on which it is block diagonal with blocks
γg
I|g|−1 , g ∈ E.
λg
Hence the matrix of interest is
−1

G(z, γ)



GE (z, γ)−1 CE,−β (z, γ) 0
C−β (z, γ) =
0
0

where the blocks of GE are
T
Πg (X T X + ǫI)−1 ΠTh Vh,zh ,
GE,gh (z, γ) = Vg,z
g

g, h ∈ E.

The eigenvalues of the matrix which show up in the determinant are the same
as the reciprocal of the eigenvalues of the matrix with blocks
ḠE,gh (z, γ) =

λg 1/2 T
λh
Vg,zg (X T X + ǫI)−1 Vh,zh ,
γg
γh

g, h ∈ E.

(104)

Finally, the matrix H(z, γ) can be ignored if we consider X to be fixed.
Appendix E: Multiple views of the data
In this section, we consider three more algorithms that choose variables by
querying the data several times.
E.1. Top K screening
A simple way of variable selection called top K screening selects the K most
correlated features Xj with response vector y. For a given K, the unrandomized
version of such a screening selects most correlated variables (j1 , j2 , . . . , jK ) with
corresponding signs (s1 , s2 , . . . , sK ), such that
s1 XjT1 y ≥ s2 XjT2 y ≥ · · · ≥ sK XjTK y ≥

max

j6∈{j1 ,...,jK }

|XjT y|.

Selective inference in the nonrandomized setting for this problem was considered
in (Lee and Taylor, 2014).
imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

67

Here, we consider a randomized version of the optimization above yielding
the k-th most correlated variable for 1 ≤ k ≤ K, with added independent
randomization variable ωk ∼ Gk ∈ Rp−k+1 and corresponding density gk is
given by
T
maximize η T (X−A
y + ωk ) − IKk (η), where y ∼ F.
k−1
η∈Rp−k+1

(105)

Here, Ak = {j1 , j2 , . . . , jk } is the active set including the k-th step, X−Ak are
the columns of X except for the ones corresponding to the current active set
Ak , ωk is a sequence of randomizations and IKk (η) is the characteristic function
of
Kk = {η ∈ Rp−k+1 : kηk1 ≤ 1},
that is,
IKk (η) =

(

0 if η ∈ Kk
∞ otherwise.

Denote the optimal solution of (105) as
(
sk if j = jk
∗
ηk,j =
0
otherwise,

(106)

where j ∈ {1, . . . , p} \ Ak−1 = Ack−1 10 and jk = argmax XjT y + ωk,j . The
j∈Ack−1

subgradient equation in k-th step leads to a reconstruction map for the randomization given by
T
ωk = φk (y, zk ) = −X−A
y + zk ,
k−1

constraining sub-differential zk ∈ Rp−k+1 to
zk ∈ ∂IKk (ηk∗ ) = {c · u : u ∈ Rp−k+1 , ujk = sk , |uj | ≤ 1 ∀j ∈ Ack−1 , c > 0}.
Conditioning on the selection event of choosing the K most correlated variables with their corresponding signs

K
Y

n
K
Rp−k+1 : sign(XjTk y + ωk ) = sk ,
Ê{(sk ,jk )}K
=
∈
R
×
y,
{ω
}
k k=1
k=1
k=1


T
sk (XjTk y + ωk ) ≥ max
|X
y
+
ω
|,
k
=
1,
.
.
.
,
K
.
k
j
c
j∈Ak−1

we sample (y, z1 , . . . , zk ) from the selective sampling density proportional to
f (y) ·

K
Y

k=1



T
y ,
gk zk − X−A
k−1

(107)

10 η ∗ ∈ Rp−k+1 is indexed by this set and all p − k + 1-dimensional vectors in this and the
k
following section will be indexed by Ack−1 .

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

68

supported on
Rn ×

K
Y

∂IKk (ηk∗ ).

k=1

For logistic regression, one can replace the T statistics above with the score
statistics as described in Remark 13.
E.2. Stagewise algorithms
Instead of fully projecting out the current variables at each step as in forward stepwise, one can use an incremental approach as in a stagewise algorithm
Tibshirani (2015). The first step of a randomized version of such an algorithm
might consist of solving the problem
maximize η T (X T (y − Xα0 ) + ω1 )
η:kηk1 ≤1

(108)

with α0 = η0∗ = 0. We update
α1 = α0 + δ · η0∗
for some learning rate δ > 0. Subsequent problems are given by
maximize η T (X T (y − Xαk ) + ωk )
η:kηk1 ≤1

with solution ηk∗ and
αk = δ ·

k−1
X

(109)

ηj∗ .

j=0

After K steps, the sampler density is thus proportional to
f (y) ·

K
Y

k=1


gk zk − X T (y − Xαk−1 ) ,

(y, z1 , . . . , zk ) ∈ Rn ×
p

with K = {η ∈ R : kηk1 ≤ 1}.

K
Y

(110)

∂IK (ηk∗ ).

k=1

E.3. Screening via thresholding randomized Post-LASSO
An alternative way to screen is through the following two-stage procedure where
we use randomized LASSO as in (41) in the first stage to select the model
(E, zE ). In the second stage, we solve an unpenalized, randomized program as
in (61) with the selected predictors XE
1
ǫ1
minimize ky − XE γk22 + kγk22 − ω1T γ, ((X, y), ω1 ) ∼ F × G1 ,
|E|
2
2
γ∈R

(111)

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

Tian Harris et al./Selective sampling

69

where G1 is a known distribution on R|E| and ω1 is independent from the randomization variable ω used in the first stage. We perform a second stage of
selection based on the output of the above convex program to threshold the
resulting coefficients of γ̂((X, y), ω1 ) resulting in a further selected model
Ẽ = {i : |γ̂i ((X, y), ω1 )| > aσ}
with their signs z̃Ẽ = sign(γ̂Ẽ ((X, y), ω1 )), where a is a constant and σ is again
the scaling which can be estimated by the noise variance of the selected model
as in Section 4.6. The canonical event of interest from the two-step procedure
becomes
B(E,zE ,Ẽ,z̃

Ẽ )

= {(βE ,u−E , γ) : diag(zE )βE > 0, ku−E k∞ ≤ 1,
diag(z̃Ẽ )γẼ > 0, |γi | ≥ aσ ∀ i ∈ Ẽ, kγ−Ẽ k∞ < aσ}.

Now the sampling density on ((X, y), βE , u−E , γ) becomes proportional to


  

zE
β
· g1 ǫ1 γ − X T (y − Xγ)
f (X, y) · g ǫ E − X T (y − XE βE ) + λ
u−E
0

where the optimization variables (βE , u−E , γ) are restricted to B(E,zE ,Ẽ,z̃ ) .
Ẽ

imsart-generic ver. 2008/08/29 file: selective_sampler.tex date: September 20, 2016

