# Finite- and Large- Sample Inference for Model and Coefficients in High-
dimensional Linear Regression with Repro Samples

###### Abstract

This paper presents a new and effective simulation-based approach to conduct
both finite- and large- sample inference for high-dimensional linear
regression models. We develop this approach under the so-called repro samples
framework, in which we conduct statistical inference by creating and studying
the behavior of artificial samples that are obtained by mimicking the sampling
mechanism of the data. We obtain confidence sets for either the true model, a
single or any collection of regression coefficients. The proposed approach
addresses two major gaps in the high-dimensional regression literature: (1)
lack of inference approaches that guarantee finite-sample performance; (2)
lack of effective approaches to address model selection uncertainty and
provide inference for the underlying true model. We provide both finite-sample
and asymptotic results to theoretically guarantee the performance of the
proposed methods. Besides enjoying theoretical advantages, our numerical
results demonstrate that the proposed methods achieve better coverage with
smaller confidence sets than the existing state-of-art approaches, such as
debiasing and bootstrap approaches. We also extend our approaches to drawing
inference on functions of the regression coefficients.

00X00,

00X00,

00X00,

high-dimensional inference,

model selection inference,

###### keywords:

[class=MSC2010]

###### keywords:

00footnotetext: Peng Wang is Associate Professor, Department of Operations,
Business Analytics, and Information Systems, University of Cincinnati,
Cincinnati, OH 45220. Email<EMAIL_ADDRESS>Min-ge Xie is Distinguished
Professors, Linjun Zhang is Assistant Professor, Department of Statistics,
Rutgers University, Piscataway, NJ 08854. Emails<EMAIL_ADDRESS>and
<EMAIL_ADDRESS>The authors acknowledge support in part from grants:
NSF-DMS DMS1812048, DMS2015373, DMS2027855, DMS2015378

##  1 Introduction

High-dimensional linear regression plays an important role in modern
statistics, with applications ranging from signal processing [52] to
econometrics [4, 36] to bioinformatics[51, 32]. There have been a number of
approaches to this problem in the past three decades, in the early part of
which researchers focused more on estimation problems such as coefficient
estimation [12, 46, 15] and support recovery [58, 64]. More recently, since
the seminal work [55, 42, 35] on the debiased LASSO, the more difficult task
of inference problems comes to the central stage, including confidence
interval construction [55, 42, 35], multiple testing of regression
coefficients [37], and post-selection inference [45, 47].

Despite a large amount of literature on inference for the high-dimensional
linear model, several important open problems remain. First, all the results
in the literature on high-dimensional inference are asymptotic, relying on the
assumption that the sample size goes to infinity, and there are no theories
concerning the performance of these procedures under the finite-sample
setting. On many occasions, especially in the high-dimensional setting with
limited sample size, the performance of these asymptotic procedures are
empirically unsatisfactory, therefore a procedure with guaranteed finite-
sample performance is desirable. Second, most existing work focus on parameter
(more specifically, regression coefficients) estimation and its inference. The
uncertainty quantification for model selection is mostly omitted in the
literature, even though it is important to account for such uncertainty for
the overall inference. This is partly due to the challenges arising from the
discrete nature of the model support. As a result, conventional inference
tools built for continuous parameters, such as the central limit theorem and
bootstrap theorems, no longer apply. Third, the post-selection inference
framework only makes statistical inference for covariates selected by some
model selection procedure. Therefore, if a covariate is significant but not
selected, no inference will be made, making some important signals being
missed. To solve these open problems, we develop a repro samples method in
this article. The proposed method quantifies both the uncertainty in model
selection and that in estimation of regression coefficients or their
functions. Although our work focuses primarily on finite-sample performances,
we also provide large-sample results when applicable.

Let us consider the high-dimensional linear regression problem in which we
observe an $n\times 1$ response vector ${\bf y}_{obs}$ with an $n\times p$
covariate matrix ${\bf X}$, where the dimension $p$ could be much larger than
the sample size $n$. Suppose that the observed ${\bf y}_{obs}$ relates only to
a subset of covariates indexed by $\tau_{0}$ with

$\displaystyle\mathbf{y}_{obs}=\mathbf{X}\bm{\beta}_{0}^{full}+\epsilon^{rel}=\mathbf{X}_{\tau_{0}}\bm{\beta}_{0}+\mathbf{\epsilon}^{rel},$

where $\epsilon^{rel}\sim N(0,\sigma_{0}^{2}\mathbf{I}_{n})$ is the
realization of the error term associated with the observed data ${\bf
y}_{obs}$, $\mathbf{I}_{n}$ is a $n\times n$ identity matrix,
$(\tau_{0},{\bm{\beta}}_{0},\sigma_{0}^{2})$ are unknown model parameters,
$\mathbf{X}=(\mathbf{X}_{\tau_{0}},\mathbf{X}_{\tau^{C}_{0}})$,
$\bm{\beta}^{full}_{0}=(\bm{\beta}^{\top}_{0},\mathbf{0}^{\top}_{\tau_{0}^{C}})^{\top}$
and $\tau_{0}^{C}=\\{1,\ldots,p\\}\setminus\tau_{0}$. We can re-express the
above model as

${\bf
y}_{obs}=\mathbf{X}\bm{\beta}_{0}^{full}+\sigma_{0}\mathbf{u}^{rel}=\mathbf{X}_{\tau_{0}}\bm{\beta}_{0}+\sigma_{0}\mathbf{u}^{rel},$
(1)

where $\mathbf{u}^{rel}$ is the realization of the standardized error term
$\mathbf{u}\sim N(0,I_{n})$ and, compared with $\epsilon^{rel}$, it is free of
the unknown $\sigma_{0}$. Corresponding to model (1) for the observed ${\bf
y}_{obs}$ and under the frequentist inference framework, there is a random
sample (or population) version of data generation model

${\bf
Y}=\mathbf{X}\bm{\beta}_{0}^{full}+\sigma_{0}\mathbf{U}=\mathbf{X}_{\tau_{0}}\bm{\beta}_{0}+\sigma_{0}\mathbf{U},\quad\hbox{$\mathbf{U}\sim
N(0,I_{n})$,}$ (2)

of which model (1) is a realization.

To make inference for the parameter $\bm{\beta}_{0}^{full}$, existing
approaches are limited by having to rely on either higher-order estimation and
asymptotics, or certain model selection procedure and its asymptotic
consistency. Accurate estimation of higher-order parameters such as variance-
covariance matrix or its inverse is a very challenging problem itself, and the
errors from the estimations might increase the width of the confidence sets
for $\bm{\beta}_{0}^{full}$. In addition, model selection consistency usually
requires strong conditions, see, e.g., [58, 6]. Also, importantly, these
methods rely on asymptotic properties without any guaranteed finite-sample
performance. In fact, we have frequently observed in our empirical studies
that the finite-sample performances of these large sample approaches are
unsatisfactory, especially when the dimension of the design matrix is high.

Moreover, inference for the true model $\tau_{0}$ has not been well studied in
high-dimensional settings, in contrast to the fact that estimation of
$\tau_{0}$, or model selection, has been deeply studied since the early
research of high-dimensional statistics. This is in large part due to the
challenges from the nature of the problem: model supports are distinctive,
unordered, and discrete objects; therefore, theories and methods developed for
continuous or ordered parameters, such as the central limit theorem and the
bootstrap method, no longer apply for inference of $\tau_{0}.$ Although there
exist works on model confidence sets that tackle the discrete model space [31,
26], the works are all in low-dimensional settings and require that the model
space is small, which apparently does not hold in the high-dimensional regime.
Moreover, the works of [31, 26] involve sequential testing procedures that do
not directly extend to the high-dimensional settings. In the end, constructing
model confidence sets for high-dimensional models remains a difficult open
problem.

This paper proposes a repro samples approach, which allows us to make
inference not only for the regression coefficients $\bm{\beta}_{0}^{full}$ but
also for the true model $\tau_{0}$. The approach directly provides us with the
desired confidence sets without having to estimate $(\tau_{0},\bm{\beta}_{0})$
or any other model parameters. The essential idea is to generate “synthetic
(fake) residuals” by sampling from the residual distribution and find
candidate models under which the synthetic residuals can generate an
artificial $\mathbf{y}$ that resembles the observed $\mathbf{y}_{obs}$. After
that, we construct both the confidence sets for $\tau_{0}$ and those for any
subset or function of $\bm{\beta}_{0}$, by checking each candidate model and
see whether it is likely to reproduce $\mathbf{y}_{obs}$. The checking
criterion is facilitated by a so called nuclear mapping function. To handle
the nuisance parameters, we propose schemes to compute the nuclear mapping
function, such that its distribution or conditional distribution is free of
the nuisance parameters.

Compared to existing approaches such as bootstrap and other related large
sample approaches, the proposed approach has several advantages: first, the
confidence sets have guaranteed finite-sample performance; the corresponding
large-sample supporting theories are also available; second, numerical studies
suggest that our proposed approach efficiently constructs tight confidence
sets, both for the true model $\tau_{0}$ and for subsets of
$\bm{\beta}^{full}_{0}$, and improved performance over the state-of-art
debiased method [35, 55] on the inference of a single coefficient is also
documented; third, our finite-sample theoretical results do not require any
conditions on the separation or strength of signals, such as the $\beta$-min
condition [58, 38], to achieve the desired coverage rate, although weaker
signals would demand a larger computational cost.

###  1.1 Contributions

To summarize, this paper has the following five major contributions.

  1. 1.

We propose the repro samples framework to effectively construct model
confidence set and quantify model selection uncertainty in the high-
dimensional linear regression model. To the best of our knowledge, this is the
only existing approach that is computationally efficient and that achieves the
goal of constructing a performance-guaranteed model confidence set without
splitting the data or predetermining candidate models.

  2. 2.

We develop a novel inference procedure for regression coefficients. Contrary
to other existing methods, our approach does not rely on the estimation of
covariance matrices or a consistent model selection procedure.

  3. 3.

Theoretically, we show that the proposed inference procedures for both the
model and the regression coefficients are able to achieve finite-sample
coverage. Although there have been many studies on high-dimensional inference,
most literature focus only on asymptotic properties. We believe the proposed
method is the first approach that guarantees the coverage for finite samples.
Additionally, our theory suggests a complementary effect between computational
power and sample size: one can achieve valid coverage as long as either is
sufficiently large.

  4. 4.

We do not need to impose the standard conditions that high-dimensional
literature typically require to obtain a consistent estimation, such as the
restricted isometry property or restricted eigenvalue conditions. Neither do
we need to require a uniform signal strength condition, which is usually
necessary for consistent model selection.

  5. 5.

Finally, through extensive numerical simulations, we show that our proposed
method is able to produce smaller confidence sets than those of the state-of-
art debiased LASSO estimators [55, 35]. In particular, due to the finite-
sample validity guarantee of our method, our proposed method is able to
achieve the desired coverage in the small sample size region, while other
methods are not.

###  1.2 Connections with other simulation-based procedures

The repro samples framework is connected to other modern simulation-based
procedures. As a technique to reproduce synthetic (“fake”) sample, it is in
principle akin to the bootstrap [22, 13]. The difference is that bootstrap
generates synthetic data from empirical distribution of estimated residuals
that are based on model selection and regression coefficients estimation,
while repro samples does so from the true residual distribution without doing
any estimation preceding the repro samples procedure. More importantly, the
bootstrap approach relies on large-sample central limit theorem (or large
deviation theorem) to justify its validity, while the repro samples can
directly provide finite-sample inference without invoking any large-sample
theory. Besides connections with the bootstrap, the repro samples method also
inherits the “matching scheme” of the approximate Bayesian computation (ABC)
approach [3], the inferential models [39] and the generalized fiducial
inference [30]. Nevertheless, unlike these Bayesian or fiducial approaches,
the proposed framework is a frequentist approach, and we do not require either
a prior, or a summary (sufficient) statistics, or an approximate match within
an ad-hoc tolerance. Finally, our framework is remotely related to the knock-
off approach [2, 11] where artificial data are used to control the false
discovery rate. See further discussions and comparisons in [53].

###  1.3 Related work

There has been much efforts in recent years developing inference procedures
for regression coefficients of $\bm{\beta}_{0}^{full}$ or functions of
$\bm{\beta}_{0}^{full}$ in high-dimensional linear regression models. On
inference of a single parameter among $\bm{\beta}_{0}^{full}$, [57, 48, 35]
propose the debiased LASSO estimator and develop its asymptotic distribution.
Another class of methods utilize Neyman-Orthogonal estimating equations with
accurate estimation of the nuisance parameters [5, 16, 17, 19]. Besides
inference for single regression coefficients, inference of linear, quadratic
and more general functionals of $\bm{\beta}_{0}^{full}$ have been studied in
[1, 7, 8, 29, 33, 34, 50, 60, 61, 28, 40]. On simultaneous inference for a
subset or all of $\bm{\beta}_{0}^{full},$ [62, 41, 49] develop high-
dimensional versions of chi-square tests, and [56, 21, 18] combine bootstrap
with the debiased approach to construct simultaneous confidence sets.
Moreover, [10, 9, 42, 59, 61] investigate constructing confidence regions
(sets) for $\bm{\beta}_{0}^{full}$. Also see [20] for a review on high-
dimensional statistical inference. The current paper provides solutions to the
same inference problems, plus inference problems on model selection. It also
provides previously unavailable finite-sample results.

A recent work by the authors [53] provided a repro samples framework for
statistical inference under a general setup in which the number of parameters
$p$ is less than the number of observations $n$. The current paper focuses on
the high-dimensional $p\gg n$ case that was not discussed in [53]. New
procedures and theoretical results with conditions tailored to high-
dimensional models that guarantee the performance of the proposed method in
both finite and large sample cases are developed.

###  1.4 Notation and identifiability condition

In this paper, for any $k\in\mathbb{N}^{+}$, we use $[k]$ to denote the set
$\\{1,2,\cdots,k\\}$. For a vector $\mathbf{v}$, let $v_{i}$ be the $i$-th
coordinate of $\mathbf{v}$. For a set $S$, we use $|S|$ to denote the
cardinality of $S$. For two positive sequences $\\{a_{k}\\}$ and
$\\{b_{k}\\}$, we write $a_{k}=O(b_{k})$ (or $a_{n}\lesssim b_{n}$), and
$a_{k}=o(b_{k})$, if $\lim_{k\rightarrow\infty}(a_{k}/b_{k})<\infty$ and
$\lim_{k\rightarrow\infty}(a_{k}/b_{k})=0$, respectively. We also use
$\mathbb{P}$ for probability and $\mathbb{E}$ for expectation in general.
Sometimes, to stress the source of randomness, we add a subscription and use
$\mathbb{P}_{Q}$ and $\mathbb{E}_{Q}$ to indicate the probability and
expectation are taken over the random variable $Q$. We use $\hat{\mathbb{P}}$
and $\hat{\mathbb{E}}$ for empirical probability and expectation. For a
$p$-dimensional vector $\bm{\beta}$ and an indices set
$\tau\subset\\{1,\ldots,p\\}$, we use $\bm{\beta}_{\tau}$ to denote the sub-
vector of $\bm{\beta}$, containing all the entries of $\bm{\beta}$ that are
associated with the indices in $\tau$. The model space ${\cal
M}=2^{\\{1,\ldots,p\\}}$ is of size $2^{p}$ and it is the collection of all
possible $\tau$ models. For a matrix $M\in\mathbb{R}^{m\times n}$, we use
$\operatorname{span}(M)$ to denote the vector space spanned by the columns of
$M$: $\operatorname{span}(M)=\\{Mv:v\in\mathbb{R}^{n}\\}$. We also call
$M(M^{\top}M)^{-1}M^{\top}$ the projection matrix of $M$.

Finally, for $\tau_{0}$ defined in (1), there might be another model
$\tilde{\tau}_{0}$ and corresponding coefficients
$\bm{\beta}_{\tilde{\tau}_{0}}$ such that
$\mathbf{X}_{\tau_{0}}\bm{\beta}_{0}=\mathbf{X}_{\tilde{\tau}_{0}}\bm{\beta}_{\tilde{\tau}_{0}}$.
Then it is not possible to tell apart $\tau_{0}$ and $\tilde{\tau}_{0}$ even
when knowing both ${\bf y}_{obs}$ and the realization $\mathbf{u}^{rel}$,
since ${\bf
y}_{obs}=\mathbf{X}_{\tau_{0}}\bm{\beta}_{0}+\sigma_{0}\mathbf{u}^{rel}=\mathbf{X}_{\tilde{\tau}_{0}}\bm{\beta}_{\tilde{\tau}_{0}}+\sigma_{0}\mathbf{u}^{rel}.$
We refer this as an identifiability issue. To address this issue and uniquely
define $\tau_{0}$, we follow the convectional practice in the high-dimensional
regression literature [6] to re-define $\tau_{0}$ as the smallest model among
the set
$\\{\tau\in\mathcal{M}|\mathbf{X}_{\tau}\bm{\beta}_{\tau}=\mathbf{X}_{\tau_{0}}\bm{\beta}_{0},\hbox{for
some }\bm{\beta}_{\tau}\\}$:

$\displaystyle\tau_{0}^{\prime}=\underset{\small\big{\\{}\tau|\mathbf{X}_{\tau}\bm{\beta}_{\tau}=\mathbf{X}_{\tau_{0}}\bm{\beta}_{0}\big{\\}}}{\rm
argmin}|\tau|.$ (3)

Throughout the paper, we assume that the true model $\tau_{0}^{\prime}$
defined in (3) is unique, which we refer as the identifiability condition. For
notational simplicity, we will still refer $\tau_{0}^{\prime}$ as $\tau_{0}$
and $\bm{\beta}_{\tau_{0}^{\prime}}$ as $\bm{\beta}_{0}$ in the remaining of
the paper.

Under the identifiability condition, we follow [44] to define the degree of
separation between model $\tau_{0}$ and other models of equal or smaller model
sizes as

$\displaystyle
C_{\min}=\min_{\\{\tau:\tau\neq\tau_{0},|\tau|\leq|\tau_{0}|\\}}\frac{1}{n\max(|\tau_{0}\setminus\tau|,1)}\|\mathbf{X}_{0}\bm{\beta}_{0}-\mathbf{X}_{\tau}\bm{\beta}_{\tau}\|_{2}^{2}.$

By the identifiability condition, $C_{\min}>0$. We remark that the notion
$C_{\min}$ is similar to $\bm{\beta}_{min}$ of the $\bm{\beta}$-min condition
in the literature. However, unlike existing literature, we do not impose any
assumption on $C_{\min}$ other than $C_{\min}>0$. We only use it as a measure
of separation.

###  1.5 Organization

The paper is organized as follows. Section 2 provides a repro sampling
framework for the inference of model selection in high-dimensional linear
models, and provides both finite-sample and large-sample guarantees. Section 3
studies the inference of regression coefficients, including
$\bm{\beta}_{0}^{full}$ and any linear transformation of
$\bm{\beta}_{0}^{full}$, and further extends to the inference of any function
of $\bm{\beta}_{0}^{full}$. Section 4 provides numerical illustrations of our
proposed method, and compare the coverage and the size of the constructed
confidence sets with the state-of-art debiased LASSO method. In Section 5, we
present technical proofs of our main results. Section 6 concludes the paper
with a discussion of our results and future research directions. The proofs of
other results and technical lemmas are deferred to the Appendix.

##  2 Inference for the true model $\tau_{0}$

In this section, we construct model confidence set for high-dimensional linear
regression models using a repro samples method. Since the model space is
discrete and of a large size $2^{p}$, to reduce the search space and improve
computational efficiency, we first introduce a notion of candidate models and
develop a data-drive method to search for the candidate models in Section 2.1.
In Section 2.2, we will then develop an inference procedure for the true model
$\tau_{0}$ with provable guarantees. We will also utilize the data-driven
candidate set in the later sections when we conduct inference for subsets of
regression coefficients, as well as inference for functions of regression
coefficients.

###  2.1 Searching for candidate models

To introduce the insight and intuition of our proposed procedure for the
search of candidate models, we first investigate how we can recover the true
model in an ideal (but unrealistic) case, that is, if we had known the
unobserved realization of the error terms $\mathbf{u}^{rel}.$ In particular,
Lemma 1 below suggests that $\tau_{0}$ defined in (3) can be expressed in
terms of the given realization $({\bf y}_{obs},\mathbf{u}^{rel})$ through an
optimization statement. We provide a proof of the following lemma in the
Appendix.

######  Lemma 1.

Let $\mathbf{H}_{\tau}$ be the projection matrix of $\mathbf{X}_{\tau}$ and
$\mathbf{H}_{\tau,\mathbf{u}^{rel}}$ be the projection matrix of
$(\mathbf{X}_{\tau},\mathbf{u}^{rel}).$ Let
$\gamma^{2}_{(\mathbf{u}^{rel},\tau_{0})}=1-\min_{\\{\tau:|\tau|<|\tau_{0}|\\}}\frac{\|(I-\mathbf{H}_{\tau,\mathbf{u}^{rel}})\mathbf{X}_{0}\bm{\beta}_{0}\|^{2}}{\|(I-\mathbf{H}_{\tau})\mathbf{X}_{0}\bm{\beta}_{0}\|^{2}}<1.$
Then given $\mathbf{u}^{rel}$, $\tau_{0}$ defined in (3) satisfies

$\displaystyle\tau_{0}$ $\displaystyle=$ $\displaystyle\underset{\tau}{\rm
argmin}\left\\{\min_{\bm{\beta}_{\tau},\sigma}\left\\{\|\mathbf{y}_{obs}-\mathbf{X}_{\tau}\bm{\beta}_{\tau}-\sigma\mathbf{u}^{rel}\|^{2}_{2}+\lambda|\tau|\right\\}\right\\},$
(4)

and

$\displaystyle(\tau_{0},\bm{\beta}_{0},\sigma_{0})$ $\displaystyle=$
$\displaystyle\underset{\tau,\bm{\beta}_{\tau},\sigma}{\rm
argmin}\left\\{\|\mathbf{y}_{obs}-\mathbf{X}_{\tau}\bm{\beta}_{\tau}-\sigma\mathbf{u}^{rel}\|^{2}_{2}+\lambda|\tau|\right\\},$

for any $0<\lambda<n[1-{\gamma^{2}_{(\mathbf{u}^{rel},\tau_{0})}}]C_{\min}.$

In practice, however, we never know the underlying realization
$\mathbf{u}^{rel}$, making it infeasible to directly apply Lemma 1.
Nonetheless, we can generate a large number of artificial samples
$\mathbf{u}^{*}\sim\mathbf{U}$, with which we construct a set of candidate
models for $\tau_{0}.$ More specifically, we generate $d$ Monte Carlo copies
$\mathbf{u}_{1}^{*},\dots,\mathbf{u}^{*}_{d}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mathbf{U}$. Then instead of solving (4) with the realized
$\mathbf{u}^{rel},$ we obtain

$\displaystyle\hat{\tau}_{b}=\underset{\tau}{\rm
argmin}\left\\{\min_{\bm{\beta}_{\tau},\sigma}\left\\{\|\mathbf{y}_{obs}-\mathbf{X}_{\tau}\bm{\beta}_{\tau}-\sigma\mathbf{u}^{*}_{b}\|^{2}_{2}+\lambda|\tau|\right\\}\right\\},$
(5)

for each $\mathbf{u}^{*}_{b},b=1,\dots,d.$ After that, we collect all
$\hat{\tau}_{b}$’s to form a candidate set for $\tau_{0}$:

$\displaystyle S^{(d)}=\left\\{\hat{\tau}_{b},b=1,\dots,d\right\\}.$ (6)

Although we require $d$ to be a large number, the candidate set $S^{(d)}$ is
often of reasonable size, since the mapping function from $\mathbf{u}^{*}_{b}$
to $\hat{\tau}_{b}$ in (5) is a many-to-one mapping and many of
$\hat{\tau}_{b}$’s obtained by (5) are identical.

The only difference between (5) and (4) is that we replace $\mathbf{u}^{rel}$
with $\mathbf{u}^{*}_{b}.$ One could imagine that if some
$\mathbf{u}^{*}_{b}\approx\mathbf{u}^{rel},$ then for such
$\mathbf{u}^{*}_{b}$’s, the event $\\{\hat{\tau}_{b}=\tau_{0}\\}$ is very
likely to happen by Lemma 1. As noted previously, the mapping function in (5)
is many-to-one, we often can correctly recover $\tau_{0}$ as long as
$\mathbf{u}_{b}^{*}$ is in a neighbourhood of $\mathbf{u}^{rel}.$ The size of
such neighbourhood depends on the separation metric $C_{\min}$ and the sample
size, yet its probability measure is always positive under the identifiability
condition. As a result, as long as $d$, the number of Monte Carlo copies, is
large enough, eventually some $\mathbf{u}_{b}^{*}$ will fall in this
neighbourhood, leading to $\hat{\tau}_{b}=\tau_{0}$ and that the candidate set
$S^{(d)}$ contains the true model $\tau_{0}$. Formally, we summarize the
results in Theorems 1-2 below, in which we prove that $\mathbb{P}_{({\cal
U}^{d},\mathbf{Y})}(\tau_{0}\notin S^{(d)})\rightarrow 0,$ in each of the
following two cases: 1) the sample size $n$ is finite and $d\to\infty$; or 2)
$d$ is finite and $n\to\infty$, respectively. Here, the probability
$\mathbb{P}_{({\cal U}^{d},\mathbf{Y})}(\cdot)$ refers to the joint
distribution of $\mathbf{U}$ and ${\cal
U}^{d}=(\mathbf{U}^{*}_{1},\dots,\mathbf{U}^{*}_{d}),$ where
$\mathbf{U}^{*}_{b}\sim\mathbf{U}$ is a Monte Carlo copy of $\mathbf{U},$
$b=1,\dots,d.$ We provide the poofs of the two theorems in Section 5.

######  Theorem 1.

Suppose $n-|\tau_{0}|>4.$ For any $\delta>0$, there exists a constant
$\gamma_{\delta}>0$ such that when
$\lambda\in\big{[}n\gamma_{\delta}^{1/2}\left\\{2+2(|\tau_{0}|+1)\frac{\log(p/2)}{n}\right\\},\allowbreak
n\gamma_{\delta}^{1/4}\frac{C_{\min}}{6}\big{]},$ the finite probability bound
that the true model is not included in the model candidates set $S^{(d)}$,
obtained by (6) with the objective function (5), is as follows,

$\displaystyle\mathbb{P}_{({\cal U}^{d},\mathbf{Y})}(\tau_{0}\notin
S^{(d)})\leq\left(1-\frac{\gamma_{\delta}^{n-1}}{n-1}\right)^{d}+\delta.$ (7)

Therefore as $d\rightarrow\infty,$ $\mathbb{P}_{({\cal
U}^{d},\mathbf{Y})}(\tau_{0}\notin S^{(d)})\rightarrow\delta,$ where
$\delta>0$ is arbitrarily small.

######  Theorem 2.

Suppose
$\frac{\lambda}{n}\in\left[6\sigma_{0}^{2}\frac{(|\tau_{0}|+1)\log(p/2)}{n}+t,0.05C_{\min}\right]$,
where $t>0$ is a positive constant. Then the finite probability bound that the
true model is not included in the confidence set $S^{(d)}$, obtained by (6)
with the objective function (5) for any finite $d$ is as follows,

$\displaystyle\mathbb{P}_{({\cal U}^{d},\mathbf{Y})}(\tau_{0}\notin S^{(d)})$
$\displaystyle\leq
6\exp\left[-\frac{n}{18\sigma_{0}^{2}}\\{0.3C_{\min}-36\frac{\log(p+1)}{n}\sigma_{0}^{2}\\}\right]+3\exp\left(-\frac{nt}{3\sigma_{0}^{2}}\right)$
$\displaystyle+\exp\left\\{-nd\left(0.23-\frac{|\tau_{0}|\log(p)+2}{n}\right)\right\\},$
(8)

Therefore $\mathbb{P}_{({\cal U}^{d},\mathbf{Y})}(\tau_{0}\notin
S^{(d)})\rightarrow 0$ for any $d$ as $n\rightarrow\infty$, if
$\frac{|\tau_{0}|\log(p)}{n}<0.23$ and
$C_{\min}>120\frac{(|\tau_{0}|+2)\log(p+1)}{n}\sigma_{0}^{2}$ when $n$ is
large enough.

The two theorems above suggest two complementary driving forces of the
coverage validity: the sample size and the computation time. In the case where
the sample is limited, Theorem 1 implies that we will have valid coverage as
long as the computation time goes to infinity; in the case where the
computational resources are limited, Theorem 2 then indicates collecting
sufficient samples will result in valid coverage guarantee.

To put it succinctly, we summarize the aforementioned procedure in Algorithm 1
below.

Algorithm 1 Search of Candidate Models

Input: Design matrix $\mathbf{X}$, response vector $\bm{y}_{obs}$, the number
of Monte Carlo copies $d$.

Output: Candidate Models $S^{(d)}$.

Step 1: Simulate a large number $d$ copies of ${\bf u}^{*}\sim\mathbf{U}\sim
N(0,\mathbf{I}_{n})$. Denote the $d$ copies by ${\bf u}^{*}_{b}$,
$b=1,\ldots,d$, respectively.

Step 2: Compute $\hat{\tau}_{b,\lambda}=\underset{\tau}{\rm
argmin}\,\underset{(\bm{\beta}_{\tau},\sigma)}{\min}\left[\lambda|\tau|+\|\mathbf{y}_{obs}-\mathbf{x}_{\tau}\bm{\beta}_{\tau}-\sigma{\bf
u}^{*}_{b}\|^{2}_{2}\right]$ for $b=1,\ldots,d$ and a grid values of
$\lambda.$ For each $b$, use certain selection criteria to pick a subset of
all values of $\lambda$, denoted as $\Lambda_{b}$.

Step 3: Construct
$S^{(d)}=\left\\{\hat{\tau}_{b,\lambda}:\lambda\in\Lambda_{b},b=1,\dots,d\right\\}$.

To apply Algorithm 1 in practice, one need to consider two components. First,
researchers commonly use some criteria to determine the value of the tuning
parameter $\lambda$. In this paper, we use the extended BIC [14] to determine
$\lambda$, due to its good empirical performance and asymptotic model
selection consistency in high-dimensional settings. The second practical
component is that solving an optimization problem with a $L_{0}$ penalty
$\lambda|\tau|$ is often computationally difficult for high-dimensional data,
and researchers typically use a surrogate to replace the $L_{0}$ penalty. In
this paper, we adopt the adaptive LASSO [63] as a surrogate for the $L_{0}$
penalty in (5) because of its simplicity and convexity. One may also use other
surrogates like the truncated LASSO penalty [43], smoothly clipped absolute
deviation penalty [23], or the minimax concave penalty [54], among others.
These surrogates may enjoy different computational and theoretical properties.
However, since the implementation of penalty functions is not the focus of
this paper, we will not dive into the topic any further.

###  2.2 Construction of Model Confidence Sets

In this subsection, we construct the confidence sets for model $\tau_{0}$ by
developing a repro samples method that is tailored to our problem.

For the ease of presenting the general repro samples method, let us first
assume that we are interested in making a joint inference about
$\bm{\theta}_{0}=(\tau_{0},\bm{\beta}_{0}^{\top},\sigma_{0})^{\top}$ and
describe how a repro samples method works in this case. The idea is that for
any potential candidate value
$\bm{\theta}=(\tau,\bm{\beta}_{\tau}^{\top},\sigma)^{\top}$, we can create an
artificial data
$\mathbf{y}^{*}=\mathbf{X}_{\tau}\bm{\beta}_{\tau}+\sigma\mathbf{u}^{*},$
where $\mathbf{u}^{*}\sim N(0,\mathbf{I}_{n}).$ This artificial data is
referred to as a copy of repro sample. If $\mathbf{u}^{*}$ is close to
$\mathbf{u}^{rel}$ and $\bm{\theta}$ is equal or close to $\bm{\theta}_{0}$,
then we expect $\mathbf{y}^{*}$ and $\mathbf{y}_{obs}$ would be equal or
close. Inversely, for a potential value $\bm{\theta}$, if we can find a
$\mathbf{u}^{*}$ likely matching $\mathbf{u}^{rel}$ such that $\mathbf{y}^{*}$
matches $\mathbf{y}_{obs}$ (i.e., $\mathbf{y}_{obs}=\mathbf{y}^{*}$), then we
keep this $\bm{\theta}$ as a potential estimate of $\bm{\theta}_{0}$.
Mathematically, we define a level-$\alpha$ confidence set for
$\bm{\theta}_{0}$ as

$\Gamma^{\theta}_{\alpha}(\mathbf{y}_{obs})=\left\\{\bm{\theta}:\exists\mathbf{u}^{*}\
\ \mbox{s.t.
}\mathbf{y}_{obs}=G(\bm{\theta},\mathbf{u}^{*}),T(\mathbf{u}^{*},\bm{\theta})\in
B_{\alpha}(\bm{\theta})\right\\}.$ (9)

Here, the function $T(\cdot,\cdot)$ is referred to as a nuclear mapping
function and the set $B_{\alpha}(\bm{\theta})$ is a fixed level-$\alpha$ Borel
set in $\mathbb{R}^{d}$ such that

${\mathbb{P}}_{\mathbf{U}}\left(T(\mathbf{U},\bm{\theta})\in
B_{\alpha}(\bm{\theta})\right)\geq\alpha.$ (10)

The interpretation of (9) is: for a potential value $\bm{\theta}$, if there
exists a $\mathbf{u}^{*}$ satisfying $T(\mathbf{u}^{*},\bm{\theta})\in
B_{\alpha}$ such that the artificial data
${\mathbf{y}}^{*}=\mathbf{X}_{\tau}\bm{\beta}_{\tau}+\sigma\mathbf{u}^{*}$
matches with $\mathbf{y}_{obs}$ (i.e., $\mathbf{y}^{*}=\mathbf{y}_{obs}$),
then we keep this $\theta$ in the set. The repro samples method uses
${\mathbb{P}}_{\mathbf{U}}\left(T(\mathbf{U},\bm{\theta})\in
B_{\alpha}(\bm{\theta})\right)$, for a given $\bm{\theta}$, as a way to
quantify the the uncertainty of $\mathbf{U}$ thus also the uncertainty of
$\mathbf{Y}$. Also, for any nuclear mapping function
$T(\mathbf{U},\bm{\theta})$, as long as we have a set
$B_{\alpha}(\bm{\theta})$ such that (10) holds, the set
$\Gamma_{\alpha}(\mathbf{z}_{obs})$ in (9) is a level-$\alpha$ confidence set.
Thus, the role of $T(\mathbf{U},\bm{\theta})$ under the repro samples frameork
is similar to that of a test statistic under the classical (Neyman-Pearson)
hypothesis testing framework. Besides, a good choice for
$T(\mathbf{U},\bm{\theta})$ is problem specific.

Now we extend the general method to make inference for only the unknown true
model $\tau_{0}$. First, in order to define the nuclear mapping
$T(\mathbf{u},\bm{\theta})$ for the inference for $\tau_{0},$ we write

$\displaystyle\mathbf{Y}_{\bm{\theta}}=\mathbf{X}_{\tau}\mathbf{\bm{\beta}}_{\tau}+\sigma\mathbf{U},$
(11)

and it is an artificial data generated from a given set of parameters
$\bm{\theta}=(\tau,\bm{\beta}_{\tau}^{\top},\sigma)^{\top}$ and $\mathbf{U}.$
Accordingly, let
$\mathbf{y}_{\bm{\theta}}=\mathbf{X}_{\tau}\mathbf{\bm{\beta}}_{\tau}+\sigma\mathbf{u},\mathbf{u}\sim\mathbf{U}$
be a realization of $\mathbf{Y}_{\bm{\theta}},$ or certain (realized) fake
data generated by (11). Then based on the artificial repro data
$(\mathbf{X},\mathbf{y}_{\bm{\theta}}),$ one can obtain an estimate of $\tau$,
denoted by $\hat{\tau}(\mathbf{y}_{\bm{\theta}}).$ In principle, one can
choose to use any estimator of reasonable performance. In this paper we use

$\displaystyle\hat{\tau}(\mathbf{y}_{\bm{\theta}})=\mathop{\rm
arg\,min\,}_{\tilde{\tau}\in{\cal
T},\bm{\beta}_{\tilde{\tau}}\in\mathbb{R}^{|\tilde{\tau}|}}\|\mathbf{y}_{\bm{\theta}}-\mathbf{X}_{\tilde{\tau}}\bm{\beta}_{\tilde{\tau}}\|^{2}\mbox{
s.t. }|\tilde{\tau}|\leq|\tau|.$

We define the nuclear mapping function for the inference of the true model
$\tau_{0}$ as
$T(\mathbf{u},\bm{\theta})=\hat{\tau}(\mathbf{y}_{\bm{\theta}})=\tilde{T}(\mathbf{y}_{\bm{\theta}},\tau).$

Handling nuisance parameters is a ubiquitously challenging issue in
statistical inference. Here, if we follow the general repro samples method, we
need to find a Borel set $B_{\alpha}(\bm{\theta})$ that satisfies (10), i.e.,
${\mathbb{P}}_{\mathbf{U}}\left(T(\mathbf{U},\bm{\theta})\in
B_{\alpha}(\bm{\theta})\right)=\mathbb{P}_{\mathbf{U}}(\tilde{T}(\mathbf{Y}_{\bm{\theta}},\tau)\in
B_{\alpha}(\bm{\theta}))\geq\alpha,$ for the nuclear mapping
$\tilde{T}(\mathbf{y}_{\bm{\theta}},\tau)$ defined above. However, the
distribution of $\tilde{T}(\mathbf{y}_{\bm{\theta}},\tau)$ involves the entire
parameter $\bm{\theta}=(\tau,\bm{\beta}_{\tau}^{\top},\sigma)^{\top}$,
including the nuisance parameters $\bm{\beta}_{\tau}$ and $\sigma^{2}$. The
task of directly obtaining $B_{\alpha}(\bm{\theta})$ is computationally
formidable, if not infeasible. To overcome the impact of the nuisance
parameters, we introduce below an effective conditional repro samples method
to handle the nuisance parameters and construct a confidence set for
$\tau_{0}$ with provable guarantees.

The idea of conditional repro samples method is to first find a quantity
$W(\mathbf{U},\bm{\theta})$, such that the conditional distribution of the
nuclear statistic $T(\mathbf{U},\bm{\theta})$ given
$W(\mathbf{U},\bm{\theta})=\mathbf{w}$ is free of the nuisance parameters
$(\bm{\beta},\sigma)$. Assume we have such $W(\mathbf{U},\bm{\theta})$; we
will discuss how to obtain $W(\mathbf{U},\bm{\theta})$ for our purpose later
in the section. Then, based on the conditional distribution of
$T(\mathbf{U},\bm{\theta})|W(\mathbf{U},\bm{\theta})$, we construct a Borel
set $B_{\alpha}(\tau,\mathbf{w})$ that depends on $\mathbf{w}$, the value of
the random quantity $\mathbf{W}(\mathbf{U},\bm{\theta})$, but not on
$(\bm{\beta},\sigma)$, such that

$\displaystyle\mathbb{P}\\{T(\mathbf{U},\bm{\theta})\in
B_{\alpha}(\tau,\mathbf{w})|W(\mathbf{U},\bm{\theta})=\mathbf{w}\\}=\mathbb{P}\\{\tilde{T}(\mathbf{Y}_{\bm{\theta}},\tau)\in
B_{\alpha}(\tau,\mathbf{w})|W(\mathbf{U},\bm{\theta})=\mathbf{w}\\}\geq\alpha.$
(12)

Accordingly, the marginal probability $\mathbb{P}(T(\mathbf{U},\bm{\theta})\in
B_{\alpha}(\tau,W(\mathbf{U},\bm{\theta}))\geq\alpha.$

Now, instead of directly following (9), we construct a subset in the model
space ${\cal M}$:

$\displaystyle\Gamma^{\tau}_{\alpha}(\mathbf{y}_{obs})=\\{\tau\in\mathcal{M}:\exists\mathbf{u}^{*}\
\mbox{and }\ (\bm{\beta}_{\tau},\sigma)\ \mbox{s.t.
}\;{\mathbf{y}}_{obs}=\mathbf{X}_{\tau}\bm{\beta}_{\tau}+\sigma\mathbf{u}^{*},$
$\displaystyle T(\mathbf{u}^{*},\bm{\theta})\in
B_{\alpha}\left(\tau,W(\mathbf{u}^{*},\bm{\theta})\right),\bm{\theta}=(\tau,\bm{\beta}_{\tau},\sigma)\\}.$
(13)

The following theorem suggests that the newly constructed set
$\Gamma^{\tau}_{\alpha}(\mathbf{y}_{obs})$ is a level-$\alpha$ confidence set
for the true model $\tau_{0}$, whose detailed proof is provided in Section 5.

######  Theorem 3.

Suppose the conditional distribution of $T(\mathbf{U},\bm{\theta})$ given
$W(\mathbf{U},\bm{\theta})=\mathbf{w}$ is free of $(\bm{\beta}_{\tau},\sigma)$
and the Borel set $B_{\alpha}(\tau,\mathbf{w})$ satisfies (12), then

$\displaystyle\mathbb{P}(\bm{\tau}_{0}\in\Gamma^{\tau}_{\alpha}(\mathbf{Y}))\geq\alpha,$

where the confidence set $\Gamma^{\tau}_{\alpha}(\mathbf{Y})$ is defined
following from (2.2).

The remaining task is to find the random quantity
$\mathbf{W}(\mathbf{U},\bm{\theta})$ and the Borel set
$B_{\alpha}(\tau,\mathbf{w})$ such that the conditional distribution of
$T(\mathbf{U},\bm{\theta})|W(\mathbf{U},\bm{\theta})\mathbf{w}$ is free of the
nuisance parameters $(\bm{\beta},\sigma)$ and the inequality (12) holds. Note
that we can rewrite (11) as
$\mathbf{Y}_{\bm{\theta}}=\mathbf{H}_{\tau}\mathbf{Y}_{\bm{\theta}}+(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{Y}_{\bm{\theta}}=\mathbf{H}_{\tau}\mathbf{Y}_{\bm{\theta}}+\sigma(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}$,
where
$\mathbf{H}_{\tau}=\mathbf{X}_{\tau}(\mathbf{X}_{\tau}^{\top}\mathbf{X}_{\tau})^{-1}\mathbf{X}_{\tau}^{\top}$
is the projection matrix of $\mathbf{X}_{\tau}$. Write
$\mathbf{A}_{\bm{\theta}}(\mathbf{U})=\mathbf{H}_{\tau}\mathbf{Y}_{\bm{\theta}}=\widetilde{\mathbf{A}}_{\bm{\theta}}({\mathbf{Y}}_{\bm{\theta}})$
and
$b_{\bm{\theta}}(\mathbf{U})=\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{Y}_{\bm{\theta}}\|=\widetilde{b}_{\bm{\theta}}(\mathbf{Y}_{\bm{\theta}})$.
It follows that

$\displaystyle\mathbf{Y}_{\bm{\theta}}$
$\displaystyle=\mathbf{A}_{\bm{\theta}}(\mathbf{U})+b_{\bm{\theta}}(\mathbf{U})\frac{(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}}{\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}\|}=\widetilde{\mathbf{A}}_{\bm{\theta}}(\mathbf{Y}_{\bm{\theta}})+\widetilde{b}_{\bm{\theta}}(\mathbf{Y}_{\bm{\theta}})\frac{(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}}{\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}\|}.$

In this equation, the “randomness” of $\mathbf{U}$ (and also
$\mathbf{Y}_{\bm{\theta}}$) are splitted into three pieces,
$\mathbf{A}_{\bm{\theta}}(\mathbf{U}),b_{\bm{\theta}}(\mathbf{U})$ and
$(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}/\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}\|$.
Under model equation (11),
$(\mathbf{A}_{\bm{\theta}}(\mathbf{U}),\,b_{\bm{\theta}}(\mathbf{U}))$
$=\left(\widetilde{\mathbf{A}}_{\bm{\theta}}({\mathbf{Y}}_{\bm{\theta}}),\,\widetilde{b}_{\bm{\theta}}(\mathbf{Y}_{\bm{\theta}})\right)$
is a sufficient statistic and the last piece
$(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}/\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}\|$
is an ancillary statistic that is free of the nuisance parameters
$(\bm{\beta}_{\tau},\sigma^{2})$. Based on this partition, we define for our
purpose
$\mathbf{W}(\mathbf{U},\bm{\theta})=\left(A_{\bm{\theta}}(\mathbf{U}),b_{\bm{\theta}}(\mathbf{U})\right)=\left(\widetilde{\mathbf{A}}_{\bm{\theta}}({\mathbf{Y}}_{\bm{\theta}}),\,\widetilde{b}_{\bm{\theta}}(\mathbf{Y}_{\bm{\theta}})\right)=\widetilde{\mathbf{W}}(\mathbf{Y}_{\bm{\theta}},\bm{\theta})$.
It follows immediately that the conditional distribution of
$\mathbf{Y}_{\bm{\theta}}\mid\mathbf{W}(\mathbf{U},\bm{\theta})=\mathbf{w}$ is
free of $(\bm{\beta}_{\tau},\sigma^{2}).$ So is the conditional probability
mass function of
$\tilde{T}(\mathbf{Y}_{\bm{\theta}},\tau)\mid\mathbf{W}(\mathbf{U},\bm{\theta})=\mathbf{w}$,

$\quad
p_{(\mathbf{w},\tau)}(\tau^{\prime})={\mathbb{P}}_{\mathbf{U}|\mathbf{w}}\left\\{\tilde{T}(\mathbf{Y}_{\bm{\theta}},\tau)=\tau^{\prime}\bigg{|}\mathbf{W}(\mathbf{U},\bm{\theta})=\mathbf{w}\right\\},\quad$
(14)

for any $\tau^{\prime}\in{\cal T}$. Note that, when given $(\mathbf{w},\tau)$,
we can use the model equation (11) to generate many copies of
$\mathbf{Y}_{\bm{\theta}}$ by repeated draws from $\mathbf{U}$. The
conditional probability mass function in (14) can be obtained by a Monte-Carlo
method.

We define the Borel set $B_{\alpha}(\tau,\mathbf{w})$ as

$\displaystyle
B_{\alpha}(\tau,\mathbf{w})=\left\\{\tau^{*}\in\mathcal{M}:\sum_{\\{\tau^{\prime}:\,p_{(\mathbf{w},\tau)}(\tau^{\prime})\leq
p_{(\mathbf{w},\tau)}(\tau^{*})\\}}p_{(\mathbf{w},\tau)}(\tau^{\prime})\geq
1-\alpha\right\\}.$ (15)

In Section 5, we prove that the conditional probability

$\displaystyle{\mathbb{P}}_{\mathbf{Y}_{\bm{\theta}}|\mathbf{w}}\left\\{\tilde{T}(\mathbf{Y}_{\theta},\tau)\in
B_{\alpha}(\tau,\mathbf{w})\big{|}\widetilde{\mathbf{W}}(\mathbf{Y}_{\bm{\theta}},\tau)=\mathbf{w}\right\\}\geq\alpha.$
(16)

It follows that, marginally,
${\mathbb{P}}_{\mathbf{Y}_{\bm{\theta}}}\left\\{\tilde{T}(\mathbf{Y}_{\theta},\tau)\in
B_{\alpha}\left(\tau,\widetilde{\mathbf{W}}(\mathbf{Y}_{\bm{\theta}},\tau)\right)\right\\}\geq\alpha.$
By (2.2) and using the candidate set $S^{(d)}$, we define

${\bar{\Gamma}}^{\tau}_{\alpha}(\mathbf{y}_{obs})=\Gamma^{\tau}_{\alpha}(\mathbf{y}_{obs})\bigcap
S^{(d)}=\left\\{\tau\in S^{(d)}:\tilde{T}(\mathbf{y}_{obs},\tau)\in
B_{\alpha}\left(\tau,\widetilde{\mathbf{W}}(\mathbf{y}_{obs},\tau)\right)\right\\}.$
(17)

To obtain the confidence set, we use a Monte-Carlo algorithm to compute the
conditional probability in (14). We summarize the procedure of constructing
the above model confidence set in Algorithm 2, with the size of the Monte-
Carlo simulations $J.$ Theorem 4 below states that
${\bar{\Gamma}}^{\tau}_{\alpha}(\mathbf{y}_{obs})$ is a level-$\alpha$
confidence set for $\tau_{0}$ with a guaranteed finite-sample coverage rate,
if we can repeatedly generate Monte Carlo samples. Theorem 5 states that even
when the number of Monte Carlo copies $d$ is limited,
${\bar{\Gamma}}^{\tau}_{\alpha}(\mathbf{y}_{obs})$ is still a level-$\alpha$
confidence set for $\tau_{0}$ if we have a large sample size $n.$

######  Theorem 4.

Under the conditions in Theorem 1, for any finite $n$ and $p$, and arbitrarily
small $\delta>0$ the coverage probability of model confidence set
${\bar{\Gamma}}^{\tau}_{\alpha}(\mathbf{y}_{obs})$ constructed above is
$\mathbb{P}_{\mathbf{Y}}\left\\{\tau_{0}\in{\bar{\Gamma}}^{\tau}_{\alpha}(\mathbf{Y})\right\\}\geq\alpha-\delta-o(e^{-c_{1}d})$
as $d\rightarrow\infty$ for some $c_{1}>0$.

######  Theorem 5.

Under the conditions in Theorem 2, for any finite $d$, the coverage
probability of model confidence set
${\bar{\Gamma}}^{\tau}_{\alpha}(\mathbf{y}_{obs})$ constructed above is
$\mathbb{P}_{\mathbf{Y}}\left\\{\tau_{0}\in{\bar{\Gamma}}^{\tau}_{\alpha}(\mathbf{Y})\right\\}\geq\alpha-o(e^{-c_{2}n})$
as $n\rightarrow\infty$ for some $c_{2}>0$.

Before the end of this subsection, we put together the following Algorithm 2
to construct the level-$\alpha$ model confidence set
${\bar{\Gamma}}^{\tau}_{\alpha}(\mathbf{y}_{obs})$ for $\tau_{0}$.

Algorithm 2 Confidence set construction for $\tau_{0}$

Input: Design matrix $\mathbf{X}$, response vector $\bm{y}_{obs}$, candidate
set $S^{(d)}$, simulation size $J$

Output: Confidence set of $\tau_{0}$

for $\tau_{b}\in S^{(d)}$ and $j\in{1,\dots,J}$ do

Step 1: Calculate
$\mathbf{w}_{obs}=(\mathbf{a}_{obs},b_{obs})=(\mathbf{H}_{\tau_{b}}\mathbf{y}_{obs},\|(\mathbf{I}-\mathbf{H}_{\tau_{b}})\mathbf{y}_{obs}\|).$

Step 2: Generate $\mathbf{u}_{j}^{*}\sim N(0,\mathbf{I}_{n})$, and compute
$\displaystyle\mathbf{y}_{j}^{*}=\mathbf{a}_{obs}+b_{obs}\frac{(\mathbf{I}-\mathbf{H}_{\tau_{b}})\mathbf{u}_{j}^{*}}{\|(\mathbf{I}-\mathbf{H}_{\tau_{b}})\mathbf{u}_{j}^{*}\|}.$
In addition, obtain the estimated model
$\hat{\tau}^{*}_{j}=\hat{\tau}(\mathbf{y}_{j}^{*})$ by
$\displaystyle\hat{\tau}(\mathbf{y}^{*}_{j})=\mathop{\rm
arg\,min\,}_{\tilde{\tau}\in{\cal T},\bm{\beta}_{\tilde{\tau}}\in
R^{|\tilde{\tau}|}}\|\mathbf{y}_{j}^{*}-X_{\tilde{\tau}}\bm{\beta}_{\tilde{\tau}}\|^{2}\mbox{
s.t. }|\tilde{\tau}|\leq|\tau_{b}|.$ (18)

end for1

Step 3: Estimate $p_{(\mathbf{w}_{obs},\tau_{b})}(\tilde{\tau})$ for all
$\tilde{\tau}\in\mathcal{T}$ by
$\displaystyle\hat{p}_{(\mathbf{w}_{obs},\tau_{b})}(\tilde{\tau})=\frac{1}{J}\sum_{j=1}^{J}{\bf
1}_{\\{\hat{\tau}_{j}^{*}=\tilde{\tau}\\}},$ where ${\bf 1}_{\\{\cdot\\}}$ is
the indicator function.

Step 3: Calculate
$\displaystyle\tilde{T}(\mathbf{y}_{obs},\tau_{b})=\mathop{\rm
arg\,min\,}_{\tilde{\tau}\in{\cal T},\bm{\beta}_{\tilde{\tau}}\in
R^{|\tilde{\tau}|}}\|\mathbf{y}_{obs}-X_{\tilde{\tau}}\bm{\beta}_{\tilde{\tau}}\|^{2}\mbox{
s.t. }|\tilde{\tau}|\leq|\tau_{b}|,$

Step 4: We then compute the estimated tail probability of
$\tilde{T}(\mathbf{y}_{obs},\tau_{b})$ as
$\displaystyle\hat{F}_{(\mathbf{w}_{obs},\tau_{b})}\left\\{\tilde{T}(\mathbf{y}_{obs},\tau_{b})\right\\}=\sum_{\\{\tilde{\tau}:p_{(\mathbf{w}_{obs},\tau_{b})}(\tilde{\tau})\leq
p_{(\mathbf{w}_{obs},\tau_{b})}(\hat{\tau}(\mathbf{y}_{obs}))\\}}\hat{p}_{(\mathbf{w}_{obs},\tau_{b})}(\tilde{\tau}).$

Step 5: We therefore obtain the level-$\alpha$ confidence set for $\tau_{0}$
$\displaystyle\bar{\Gamma}^{\tau}_{\alpha}(\mathbf{y}_{obs})=\left\\{\tau_{b}\in
S^{(d)}:\hat{F}_{(\mathbf{w}_{obs},\tau_{b})}\left\\{\tilde{T}(\mathbf{y}_{obs},\tau_{b})\right\\}\geq
1-\alpha\right\\}.$

##  3 Inference for regression coefficients accounting for model selection
uncertainty

We first propose a confidence set for any subset of
$\bm{\beta}_{0}^{full}=(\beta_{0,1},\dots,\beta_{0,p})^{\top}$ while
accounting for model selection uncertainty in Section 3.1. The work can be
extended to make inference for any linear transformation of
$\bm{\beta}^{full}_{0}$. We then discuss two interesting special cases of
practical importance in Section 3.2: (a) joint inference for all regression
coefficients $\bm{\beta}_{0}^{full}$; and (b) inference for single regression
coefficient $\beta_{0,i},i=1,\dots,p.$ Note that, most existing methods focus
only on one of the two special cases, and there are few effective approaches
on making inference for any subset or linear transformation of
$\bm{\beta}^{full}_{0}$ in the literature. Moreover, our work guarantees both
the finite-sample and the large-sample coverage, while the existing methods
provide only asymptotic inferences. In addition, our confidence set is a union
of multiple smaller sets. In contrast, confidence sets produced by existing
methods are single intervals or ellipsoid sets; see Section 3.2 for further
comparisons.

###  3.1 Inference for any subset of regression coefficients

Let $\bm{\beta}_{0,\Lambda},\Lambda\subset\\{1,\dots,p\\}$ be the collection
of $\beta_{0,i}$’s that are of interests, leaving the remaining parameters
$\beta_{0,i},i\not\in\Lambda$, $\sigma_{0}$ and $\tau_{0}$ as nuisance
parameters. Here our strategy is to first remove the influence of
$\beta_{0,i},i\not\in\Lambda$ and $\sigma_{0}$ by defining a nuclear mapping
function that only involves $\bm{\eta}_{\Lambda}=(\bm{\beta}_{\Lambda},\tau)$
and construct a joint confidence set for $\mathbf{\eta}_{\Lambda}$. We then
utilize the work in the previous section to handle the impact of $\tau$.

For a given $\bm{\eta}_{\Lambda}$, we define the nuclear mapping as follows

$\displaystyle
T(\mathbf{u},\bm{\eta}_{\Lambda})=\begin{cases}\frac{\mathbf{u}^{\top}\mathbf{O}_{\tau,\Lambda}\mathbf{u}}{\mathbf{u}^{\top}(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{u}/(n-|\tau|)}&\text{if
$\Lambda\cap\tau\neq\emptyset$ \& $\beta_{i}=0$ for any
$i\in\Lambda\setminus\tau$}\\\ \infty&\text{if $\beta_{i}\neq 0$ for any
$i\in\Lambda\setminus\tau$}\\\ 0&\text{if $\Lambda\cap\tau=\emptyset$ \&
$\beta_{i}=0$ for any $i\in\Lambda$ }\end{cases},$ (19)

where $\mathbf{O}_{\tau,\Lambda}$ is the projection matrix of
$(\mathbf{I}-\mathbf{H}_{\tau\setminus\Lambda})\mathbf{X}_{\Lambda\cap\tau},$
and $\mathbf{H}_{\tau\setminus\Lambda}$ is the projection matrix of
$\mathbf{X}_{\tau\setminus\Lambda}.$ We can rewrite the above nuclear mapping
as a function of
$\mathbf{y}_{\bm{\theta}}=\mathbf{X}_{\tau}\bm{\beta}_{\tau}+\sigma\mathbf{u}$:

$\displaystyle{\small\tilde{T}(\mathbf{y}_{\bm{\theta}},\bm{\eta}_{\Lambda})=\begin{cases}\frac{(\mathbf{y}_{\bm{\theta}}-\mathbf{X}_{\Lambda}\bm{\beta}_{\Lambda})^{\top}\mathbf{O}_{\tau,\Lambda}(\mathbf{y}_{\bm{\theta}}-\mathbf{X}_{\Lambda}\bm{\beta}_{\Lambda})}{(\mathbf{y}_{\bm{\theta}}-\mathbf{X}_{\Lambda}\bm{\beta}_{\Lambda})^{\top}(\mathbf{I}-\mathbf{H}_{\tau})(\mathbf{y}_{\bm{\theta}}-\mathbf{X}_{\Lambda}\bm{\beta}_{\Lambda})/(n-|\tau|)}&\text{if
$\Lambda\cap\tau\neq\emptyset$ \& $\beta_{i}=0$ for any
$i\in\Lambda\setminus\tau$}\\\ \infty&\text{if $\beta_{i}\neq 0$ for any
$i\in\Lambda\setminus\tau$}\\\ 0&\text{if $\Lambda\cap\tau=\emptyset$ \&
$\beta_{i}=0$ for any $i\in\Lambda$ }\end{cases}.}$ (20)

Since when $\Lambda\cap\tau\neq\emptyset$ and $\beta_{i}=0$ for any
$i\in\Lambda\setminus\tau,$ the distribution of the nuclear mapping defined
above is
$T(\mathbf{U},\bm{\eta}_{\Lambda})=\tilde{T}(\mathbf{Y}_{\bm{\theta}},\bm{\eta}_{\Lambda})\sim
F_{|\tau\cap\Lambda|,n-|\tau\cap\Lambda|},$ we let the Borel set be
$B_{\bm{\eta}_{\Lambda}}(\alpha)=[0,F_{|\tau\cap\Lambda|,n-|\tau\cap\Lambda|}^{-1}(\alpha)],$
such that
$\mathbb{P}(\tilde{T}(\mathbf{Y}_{\bm{\theta}},\bm{\eta}_{\Lambda})\in
B_{\bm{\eta}_{\Lambda}}(\alpha))\geq\alpha.$ Then a valid level-$\alpha$ repro
samples confidence set for
$\bm{\eta}_{0,\Lambda}=(\bm{\beta}_{0,\Lambda},\tau_{0})$ is

$\displaystyle\Gamma^{\bm{\eta}_{\Lambda}}_{\alpha}(y_{obs})=\left\\{\bm{\eta}_{\Lambda}:\tilde{T}(\mathbf{y}_{obs},\bm{\eta}_{\Lambda})\in
B_{\alpha}({\bm{\eta}_{\Lambda}})\right\\}.$

Now with both $\beta_{0,i},i\not\in\Lambda$ and $\sigma_{0}$ out of the
picture, we need to deal with the only remaining nuisance parameter,
$\tau_{0}.$ To handle the impact of $\tau_{0},$ we utilize the model candidate
set $S^{(d)}$ constructed in Section 2.1 and take a union approach. That is,
for certain $\beta_{\Lambda}$, if $(\beta_{\Lambda},\tau)$ is defined above
for any $\tau$ in the candidate set $S^{(d)},$ we then retain the
$\beta_{\Lambda}$ in the confidence set for $\beta_{0,\Lambda}.$ Specifically,

$\displaystyle\Gamma^{\bm{\beta}_{\Lambda}}_{\alpha}(y_{obs})$
$\displaystyle=\left\\{\bm{\beta}_{\Lambda}:\tilde{T}(\mathbf{y}_{obs},\bm{\eta}_{\Lambda})\in
B_{\alpha}({\bm{\eta}_{\Lambda}}),\bm{\eta}_{\Lambda}=(\bm{\beta}_{\Lambda},\tau)\text{
for some }\tau\in S^{(d)}\right\\}$ $\displaystyle=\bigcup_{\tau\in
S^{(d)}}\left\\{\bm{\beta}_{\Lambda}:\tilde{T}(\mathbf{y}_{obs},\bm{\eta}_{\Lambda})\in
B_{\alpha}({\bm{\eta}_{\Lambda}}),\bm{\eta}_{\Lambda}=(\bm{\beta}_{\Lambda},\tau)\right\\}.$
(21)

An illustration of such a confidence set is provided in Section 3.2 for the
special cases that $\Lambda=\\{1,\ldots,p\\}$; see Figure 2.

We observe that inside the union in (3.1), each set is a confidence set based
on certain model $\tau$ in the candidate set $S^{(d)}.$ Although we do not
know the true underlying model $\tau_{0}$, with Algorithm 1, we are able to
construct a candidate set of reasonable size that would include $\tau_{0}$
with high probability. This enables us to guarantee the coverage rate, both in
finite samples and asymptotically, as indicated in the following theorems. We
provide detailed proofs of the two theorems in the appendix.

######  Theorem 6.

Under the conditions in Theorem 1, for any finite $n$ and $p$ and an
arbitrarily small $\delta>0,$ the coverage probability of the confidence
interval $\Gamma^{\bm{\beta}_{\Lambda}}_{\alpha}(\mathbf{Y})$ defined in (3.1)
is
$\mathbb{P}_{\mathbf{U}}\\{\bm{\beta}_{0,\Lambda}\in\Gamma^{\bm{\beta}_{\Lambda}}_{\alpha}(\mathbf{Y})\\}\geq\alpha-\delta-o(e^{-c_{1}d})$
for some $c_{1}>0$.

######  Theorem 7.

Under the conditions in Theorem 2, for any finite $d$,
$\Gamma^{\bm{\beta}_{\Lambda}}_{\alpha}(\mathbf{Y})$ defined in (3.1) is
$\mathbb{P}_{\mathbf{U}}\\{\bm{\beta}_{0,\Lambda}\in\Gamma^{\bm{\beta}_{\Lambda}}_{\alpha}(\mathbf{Y})\\}\geq\alpha-o(e^{-c_{2}n})$
for some $c_{2}>0$.

######  Remark 1 (Extension to make inference for any linear transformation of
$\bm{\beta}^{full}_{0}$).

Let $\mathbf{L}\bm{\beta}^{full}_{0}$ be a linear transformation of
$\bm{\beta}^{full}_{0},$ where $\mathbf{L}$ is a $l\times p$ transformation
matrix. We then let

$\displaystyle\widetilde{\mathbf{L}}=\left[\begin{array}[]{c}\mathbf{L}\\\
\begin{array}[]{cc}\mathbf{I}_{(p-l)\times(p-l)},&\mathbf{0}_{(p-l)\times
l}\end{array}\end{array}\right],$

$\tilde{\bm{\beta}}_{0}^{full}=\widetilde{\mathbf{L}}\bm{\beta}_{0}^{full},$
and $\widetilde{\mathbf{X}}=\mathbf{X}\widetilde{\mathbf{L}}^{-1}.$ The
inference of $\mathbf{L}\bm{\beta}^{full}_{0}$ based on the data
$(\mathbf{y}_{obs},\mathbf{X})$ is now equivalently transformed to the
inference of a subset of $\tilde{\bm{\beta}}_{0}^{full}$ based on the
transformed data $(\mathbf{y}_{obs},\widetilde{\mathbf{X}}).$ Therefore we are
able to construct the confidence set for $\mathbf{L}\bm{\beta}^{full}_{0}$ by
applying (3.1) on $(\mathbf{y}_{obs},\widetilde{\mathbf{X}}).$ Note that one
should also derive the candidate set $S^{(d)}$ from the transformed data
$(\mathbf{y}_{obs},\widetilde{\mathbf{X}})$ using Algorithm 1.

######  Remark 2 (Use a model confidence confidence set
${\bar{\Gamma}}_{\alpha}^{\tau}(\mathbf{y}_{obs})$ to construct the confidence
set $\Gamma^{\bm{\beta}_{\Lambda}}_{\alpha}(y_{obs})$ for coefficients).

In the construction of the confidence set
$\Gamma^{\bm{\beta}_{\Lambda}}_{\alpha}(y_{obs})$ for $\bm{\beta}_{0,\Lambda}$
in (3.1), we let $\tau\in S^{(d)}$. Alternatively, we can limit $\tau$ in a
level-$\alpha_{1}$ model confidence set
$\Gamma_{\alpha_{1}}^{\tau}(\mathbf{y}_{obs})$ obtained in Section 2 using
(17). Here, $\alpha_{1}\in(\frac{1}{2},1)$ and close to $1$. Similarly, take
another $\alpha_{2}\in(\frac{1}{2},1)$, and let
$\alpha=\alpha_{1}+\alpha_{2}-1$. A modified of (3.1) is
${\Gamma}^{\bm{\beta}_{\Lambda},mod}_{\alpha}(y_{obs})=\bigcup_{\tau\in{\bar{\Gamma}}_{\alpha_{1}}^{\tau}(\mathbf{y}_{obs})}\big{\\{}\bm{\beta}_{\Lambda}:\tilde{T}(\mathbf{y}_{obs},\bm{\eta}_{\Lambda})\in
B_{\alpha_{2}}({\bm{\eta}_{\Lambda}}),\bm{\eta}_{\Lambda}=(\bm{\beta}_{\Lambda},\tau)\big{\\}}$.
We prove in Section 4 this
${\Gamma}^{\bm{\beta}_{\Lambda},mod}_{\alpha}(y_{obs})$ set is a
level-$\alpha$ confidence set for $\bm{\beta}_{0,\Lambda}$. If, for instance,
we take $\alpha_{1}=\alpha_{2}=.975$, then
$\alpha=\alpha_{1}+\alpha_{2}-1=.95$ and the above
$\Gamma^{\bm{\beta}_{\Lambda},mod}_{\alpha}(y_{obs})$ has at least $95\%$
guaranteed coverage. This construction also applies to the two special cases
discussed in the next subsection.

###  3.2 Two special examples of interest

As stated in [55, 35], we are often interested in marginal inference for a
single regression coefficient $\beta_{0,i}$ in practice. Another interesting
inference problem that has been studied in the literature is to jointly infer
all regression coefficients $\bm{\beta}_{0}^{full}$ [56, 21]. In this
subsection, we consider these interesting special cases and provide a new
solution using the repro samples method.

Inference for a single regression coefficient To obtain the repro samples
confidence set for $\beta_{0,i},$ we simplify the nuclear mapping function
defined in (19) and (20) by making $\Lambda=\\{i\\}$. For a given
$\bm{\eta}_{i}=(\beta_{i},\tau)$, it is

$\displaystyle T(\mathbf{u},\bm{\eta}_{i})$
$\displaystyle=\begin{cases}\frac{\mathbf{u}^{\top}\mathbf{O}_{\tau,i}\mathbf{u}}{\mathbf{u}^{\top}(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{u}/(n-|\tau|)}&\text{if
$i\in\tau$}\\\ \infty&\text{if $i\notin\tau$ and $\beta_{i}\neq 0$}\\\
0&\text{if $i\notin\tau$ and $\beta_{i}=0$ }\end{cases}$
$\displaystyle=\begin{cases}\frac{(\mathbf{y}_{\bm{\theta}}-\mathbf{X}_{i}\beta_{i})^{\top}\mathbf{O}_{\tau,i}(\mathbf{y}_{\bm{\theta}}-\mathbf{X}_{i}\beta_{i})}{(\mathbf{y}_{\bm{\theta}}-\mathbf{X}_{i}\beta_{i})^{\top}(\mathbf{I}-\mathbf{H}_{\tau})(\mathbf{y}_{\bm{\theta}}-\mathbf{X}_{i}\beta_{i})/(n-|\tau|)}&\text{if
$i\in\tau$}\\\ \infty&\text{if $i\notin\tau$ and $\beta_{i}\neq 0$}\\\
0&\text{if $i\notin\tau$ and $\beta_{i}=0$
}\end{cases}:=\tilde{T}(\mathbf{y}_{\bm{\theta}},\bm{\eta}_{i}),$

where $\mathbf{O}_{\tau,i}$ is the projection matrix of
$(\mathbf{I}-\mathbf{H}_{\tau_{-i}})\mathbf{X}_{i}$, and
$\tau_{-i}=\tau\setminus\\{i\\}.$ Note that for $i\in\tau,$ the nuclear
statistics $\tilde{T}(\mathbf{y}_{\bm{\theta}},\bm{\eta}_{i})$ is equivalent
to the square of t-statistics for testing $H_{0}:\beta_{0,i}=\beta_{i}$,
obtained by fitting a linear regression of $\mathbf{y}_{\bm{\theta}}$ on
$\mathbf{X}_{\tau}.$

Then following from (3.1), we obtain the confidence set for $\beta_{0,i}$

$\displaystyle\Gamma^{\beta_{i}}_{\alpha}(\mathbf{y}_{obs})$
$\displaystyle=\left\\{\beta_{i}:\tilde{T}(\mathbf{y}_{obs},\bm{\eta}_{i})\in
B_{\alpha}({\bm{\eta}_{i}}),\bm{\eta}_{i}=(\beta_{i},\tau)\text{ for some
}\tau\in S^{(d)}\right\\}$ $\displaystyle=\bigcup_{\tau\in
S^{(d)}}\left\\{\beta_{i}:\tilde{T}(\mathbf{y}_{obs},(\beta_{i},\tau))\leq
F_{1,n-|\tau|}^{-1}(\alpha)\right\\},$

where we let
$B_{\alpha}({\bm{\eta}_{i}})=B_{\alpha}({\tau})=[0,F_{1,n-|\tau|}^{-1}(\alpha)].$
Following Theorems 6 and 7, the set $\Gamma^{\beta_{i}}_{\alpha}(y_{obs})$ is
a level-$\alpha$ confidence set for $\beta_{0,i}$.

######  Remark 3 (Comparison with the debiased method).

We discuss the difference between our method and the debiased LASSO. First of
all, our proposed method offers the finite-sample coverage guarantee, while
the debiased LASSO method can only produce the asymptotic coverage rate. To
put it more clearly, the debiased LASSO method needs the sample size
$n\to\infty$ to make sure the bias, which comes from the regularized
estimation and is of order $O(|\tau_{0}|\sqrt{\log p/n})$, goes to 0. In
contrast, our method bypasses the estimation step and constructs the
confidence sets directly via the repro sampling framework, and is therefore
unbiased in nature. Second, the debiased LASSO method [35] is designed to make
inference for an individual regression coefficient. The idea of debiased LASSO
method was later generalized to make inference of functional of the
coefficients, such as co-heritability[29] and group inference statistics [28].
However, such a generalization relies on specific forms of the functional and
does not allow arbitrary functional. Our method, however, as we will show in
Remark 5, can be used to construct the confidence sets for arbitrary
functionals. Third, as we will show in the next section, we found that when
the sample size is small, the debiased LASSO method is unable to achieve
desirable coverage for large regression coefficients, while our method is able
to guarantee the coverage for both large and zero regression coefficients in
this small sample setting. We refer the readers to Section 4.2 for details.

Joint inference for all regression coefficients. Let $\Lambda=\\{1,\dots,p\\}$
and we make a joint inference all regression coefficients
$\bm{\beta}_{0}^{full}$ here. Note that $\tau\subset\Lambda$, so
$\tau\setminus\Lambda=\emptyset.$ Following from (19) and (20), the nuclear
mapping function for
$\bm{\eta}=(\bm{\beta},\tau)=\left(\left(\bm{\beta}_{\tau},\bm{0}_{\tau^{C}}\right),\tau\right)$
is

$\displaystyle T(\mathbf{u},\bm{\eta})$
$\displaystyle=\frac{\mathbf{u}^{\top}\mathbf{H}_{\tau}\mathbf{u}/|\tau|}{\mathbf{u}^{\top}(I-\mathbf{H}_{\tau})\mathbf{u}/(n-|\tau|)}$
$\displaystyle=\frac{{(\mathbf{y}_{\bm{\theta}}-\mathbf{X}_{\tau}\bm{\beta}_{\tau})}^{\top}\mathbf{H}_{\tau}(\mathbf{y}_{\bm{\theta}}-\mathbf{X}_{\tau}\bm{\beta}_{\tau})/|\tau|}{{(\mathbf{y}_{\bm{\theta}}-\mathbf{X}_{\tau}\bm{\beta}_{\tau})}^{\top}(I-\mathbf{H}_{\tau})(\mathbf{y}_{\bm{\theta}}-\mathbf{X}_{\tau}\bm{\beta}_{\tau})/(n-|\tau|)}:=\tilde{T}(\mathbf{y}_{\bm{\theta}},\bm{\eta}).$
(22)

We then let the Borel set be e
$B_{\alpha}({\bm{\eta}})=B_{\alpha}(\tau)=[0,F^{-1}_{|\tau|,n-|\tau|}(\alpha)],$
and construct the joint confidence set for $\bm{\beta}_{0}^{full}$ following
from (3.1),

$\displaystyle\Gamma_{\alpha}^{\bm{\beta}}(\mathbf{y}_{obs})$
$\displaystyle=\left\\{\bm{\beta}:\tilde{T}\left(\mathbf{y}_{obs},\bm{\eta}\right)\leq
F^{-1}_{|\tau|,n-|\tau|}(\alpha),\bm{\beta}=\left(\bm{\beta}_{\tau},\bm{0}_{\tau^{C}}\right),\bm{\eta}=(\bm{\beta},\tau)\text{
for some }\tau\in S^{(d)}\right\\}$ $\displaystyle=\bigcup_{\tau\in
S^{(d)}}\left\\{\bm{\beta}:\tilde{T}\left(\mathbf{y}_{obs},\bm{\eta}\right)\leq
F^{-1}_{|\tau|,n-|\tau|}(\alpha),\bm{\beta}=\left(\bm{\beta}_{\tau},\bm{0}_{\tau^{C}}\right),\bm{\eta}=(\bm{\beta},\tau)\right\\}.$
(23)

Again, following Theorems 6 and 7, the set
$\Gamma_{\alpha}^{\bm{\beta}}(\mathbf{y}_{obs})$ is a level-$\alpha$
confidence set for the entire vector $\bm{\beta}_{0}^{full}$.

######  Remark 4 (Visualization of the confidence set
$\Gamma_{\alpha}^{\bm{\beta}}(\mathbf{y}_{obs})$ in (3.2)).

We now use a 3-dimensional graph to present a visualization of the joint
confidence set $\Gamma_{\alpha}^{\bm{\beta}}(\mathbf{y}_{obs})$ for
$\bm{\beta}_{0}^{full}$. To do so, we consider a particular example of $p=8$
with the true model $\tau_{0}=\\{1,2\\}$, for which our candidate set contains
only three models $S^{(d)}=\big{\\{}\\{1,2,3\\},\\{1,2\\},\\{1,3\\}\big{\\}}$,
each having three or less covariates. Different from the confidence set
obtained in [56, 21], which,in this example, would typically be a
$8$-dimensional shallow disc, our confidence set
$\Gamma_{\alpha}^{\bm{\beta}}(\mathbf{y}_{obs})$ is a union of three sets, one
3-dimensional ellipsoid and two 2-dimensional ellipsoids, corresponding to
models $\\{1,2,3\\}$, $\\{1,2\\}$ and $\\{1,3\\}$, respectively. Plotted in
Figure 2 are two components: (a) a confidence curve plot of the models plotted
on the candidate model space $S^{(d)}$; and (b) the corresponding confidence
regions of the coefficients in the three candidate models. The $y$-axis of
plot (a) is the associated confidence level of each model computed via the
conditional probability in (14) , therefore the plot demonstrate the
uncertainty of the models. The figure on the right shows the level-$95\%$
confidence sets of $\bm{\beta}$ (the two blue ones) for each of the three
models in the candidate set $S^{(d)}$. It demonstrates that our algorithm
produces a union of three sets of different dimensions in this example.

Figure 1: Plot for Displaying the Joint Confidence Sets

Figure 2: (a) Confidence curve plot on $S^{(d)}$; (b) Corresponding confidence
sets (one 3-dimensional ellipsoid and two 2-dimensional ellipsoids) of the
three models in candidate set
$S^{(d)}=\big{\\{}\\{1,2,3\\},\\{1,2\\},\\{1,3\\}\big{\\}}$. In (a), the red
line instantiate the case where we aim to construct a level-$0.95$
($\alpha=0.95$) model confidence set. In this case, our $95\%$ model
confidence set for the true $\tau_{0}$ contains two models; i.e.,
${\Gamma}^{\tau}_{\alpha}(\mathbf{y}_{obs})=\big{\\{}\\{1,2,3\\},\\{1,2\\}\big{\\}}$.
In (b), a $95\%$ joint confidence set for $\bm{\beta}_{0}^{full}$ is the union
of these three confidence sets, one 3-dimensional on the
$(\beta_{1},\beta_{2},\beta_{3})$ space and two 2-dimensional ellipsoids on
the $(\beta_{1},\beta_{2})$ and $(\beta_{1},\beta_{3})$ space, respectively
(in each of the cases the remaining $\beta_{j}$’s are $0$). Also, from Remark
2, the union of the 3-dimensional ellipsoid on the
$(\beta_{1},\beta_{2},\beta_{3})$ space and the 2-dimensional ellipsoid on the
$(\beta_{1},\beta_{2})$ space (with the remaining $\beta_{j}$’s equal to $0$)
forms a $90\%$ joint confidence set for $\bm{\beta}_{0}^{full}$.

######  Remark 5 (Extension to make inference for any function of
$\bm{\beta}^{full}_{0}$).

We can extend the joint confidence set in (3.2) to obtain a repro samples
confidence set for any function of $\bm{\beta}^{full}_{0},$ say
$h(\bm{\beta}^{full}_{0})$. To put it more clearly, for each
$(\tau,\bm{\beta}_{\tau})$ in the confidence set (3.2), we collect the
corresponding function value $h(\bm{\beta}_{\tau}^{full}),$ where
$\bm{\beta}_{\tau}^{full}=(\bm{\beta}_{\tau},\bm{0}_{\tau^{C}})$ to form the
confidence set for $h(\bm{\beta}^{full}_{0}),$ i.e.
$\Gamma^{h}_{\alpha}(\mathbf{y}_{obs})=\big{\\{}h\left((\bm{\beta}_{\tau},\bm{0}_{\tau^{C}})\right):\tilde{T}\left(\mathbf{y}_{obs},\bm{\eta}\right)\leq
F^{-1}_{|\tau|,n-|\tau|}(\alpha),\bm{\eta}=(\tau,\bm{\beta}_{\tau})$ for
$\tau\in S^{(d)}\big{\\}}.$

##  4 Simulation Studies

In this section, we provide three comprehensive simulation studies to evaluate
the numerical performances of the proposed repro samples methods. Section 4.1
focuses on making inference for the unknown underlying true model discussed in
Section 2, and Section 4.2 addresses the two inference problems for regression
coefficients discussed in Section 3.2.

###  4.1 Model candidates and inference of true model $\tau_{0}$

In this subsection, we study how the proposed repro samples method performs
when we use it to obtain a data-driven set of candidate models and derive a
level-$\alpha$ confidence set for unknown true model $\tau_{0}$. We carry out
our studies on simulated data from the following two models:

  * (M1)

(Extreme High Dimension) Let Let $\bm{\beta}_{0}^{full}=(3,2,1.5,0,\dots,0)$.
The correlation between $x_{j_{1}}$ and $x_{j_{2}}$ is set to
$0.5^{|j_{1}-j_{2}|}$. We set $n=50,p=1000$ and $\sigma=1$.

  * (M2)

(Decaying Signal) Let $\bm{\beta}_{0}^{full}=(2,1.5,1,0.8,0.6,0,\dots,0).$ The
correlation between $x_{j_{1}}$ and $x_{j_{2}}$ is set to
$0.1^{|j_{1}-j_{2}|}$. We let $n=80,p=150$ and $\sigma=1$.

  * (M3)

(High Dimension, Decaying Signal) Let
$\bm{\beta}_{0}^{full}=(3,2,1.5,1,0.8,0.6,0,\dots,0).$ The correlation between
$x_{j_{1}}$ and $x_{j_{2}}$ is set to $0.1^{|j_{1}-j_{2}|}$. We let
$n=100,p=500$, and $\sigma=1$.

The first model (M1) represents an extreme high-dimensional setting with $p\gg
n.$ The second model (M2) represents a challenging case where the signals are
decaying with the weakest signal being just 0.5. We set Model (M3) by
increasing the dimension of Model (M2) to investigate the performance of the
proposed approach when both high-dimensional design matrix and weak signals
are present in the data. We also add a strong signal of $\beta_{1,0}=3$ to
expand the range of signal strengths. To evaluate the frequentist performance
of the proposed methods, we replicate the simulation for 200 times for each
model.

When applying Algorithm 1 to obtain the model candidate set, we use the
following extended BIC to choose the values of the tuning parameter in (5),

$\displaystyle\mbox{EBIC}_{b,\zeta}(\lambda)=n\log\left[\left\|\mathbf{y}_{obs}-X_{\hat{\tau}_{b,\lambda}}\hat{\beta}_{\hat{\tau}_{b,\lambda}}\right\|^{2}/n\right]+\left|\hat{\tau}_{b,\lambda}\right|\log(n)+2\zeta\log\genfrac{(}{)}{0.0pt}{0}{p}{\left|\hat{\tau}_{b,\lambda}\right|}.$

Here, $\hat{\tau}_{b,\lambda}$ is the solution to (5) with the tuning
parameter $\lambda,$ $\hat{\beta}_{\hat{\tau}_{b,\lambda}}$ is an estimation
of $\beta_{\hat{\tau}_{b,\lambda}}$ and $0\leq\zeta\leq 1$ can range between 0
and 1. To increase the efficiency of candidate models search, we pick multiple
models for each $\mathbf{u}^{*}_{b}.$ Specifically, we pick all $\lambda$’s
between $\lambda^{0}_{b}$ and $\lambda^{1}_{b}$ i.e.
$\Lambda_{b}=[\lambda^{0}_{b},\lambda^{1}_{b}],$ where
$\lambda^{\zeta}_{b}=\mathop{\rm
arg\,min\,}_{\lambda}\mbox{EBIC}_{b,\zeta}(\lambda).$ This is equivalent to
using all $0\leq\zeta\leq 1$, because $\lambda^{\zeta}_{b}$ is larger for
larger $\zeta$, and [14] showed that the model selection consistency of
extended BIC holds for some $0\leq\zeta\leq 1.$ When applying Algorithm 2 to
obtain the model confidence set, we calculate the
$\hat{\tau}(\mathbf{y}_{j}^{*})$ in (18) by obtaining the largest estimated
model that is not larger than $|\tau_{b}|$ in the adaptive lasso solution
path. Also, in our analysis, we set the number of repro samples for the
candidate set in Algorithm 1 to be $d=1,000$ for model (M1). For Model
(M2)-(M3) with weak signals, identifying the true model is a notoriously
challenging problem. In this case, we set the number of repro samples to be a
large $d=10,000$ for (M2) and $d=100,000$ for (M3). For the number of repro
samples in Algorithm 2 for calculating the distribution of the nuclear
statistics, we set $J=200$ for all the three models.

We compare the proposed repro samples approach with the residual bootstrap
approach in the current literature [27, 13], where we used three different
tuning criteria, AIC, BIC and cross-validation (CV) to select models. The
numbers of bootstrap samples are $1,000,10,000$, and $100,000$ for (M1), (M2)
and (M3) respectively, the same as the numbers of repro samples for searching
for the candidate models in our method. The bootstrap model “confidence” sets
are obtained by removing the least frequent model estimations from the
bootstrap data, with the total (cumulative) frequency of the removed models
not larger than 5%. We note that the bootstrap method here is a naive method
used in practice. Due to the discreteness of the model space and estimated
models, there is no theoretical support to the “confidence” claim that that
such a bootstrap method can get a valid level-$95\%$ model confidence set for
the true model $\tau_{0}$.

Table 1: Comparison of Performance of the Model Confidence Sets | Medthod |  Cardinality of $S^{(d)}$ |  Inclusion of $\Gamma^{\tau}_{0.95}$
---|---|---|---
Model M1: $n=50,p=1000$ | Repro samples | 2.605 (0.191) | 1.000 (0.000)
Bootstrap AIC | 215.425 (10.855) | 1.000 (0.000)
Bootstrap BIC | 146.100 ( 7.423) | 1.000 (0.000)
Bootstrap CV | 259.535 (11.891) | 1.000 (0.000)
Model M2: $n=80,p=150$ | Repro samples | 29.455 (3.080) | 0.980 (0.010)
Bootstrap AIC | 4350.850 (134.000) | 1.000 (0.000)
Bootstrap BIC | 2303.190 ( 75.708) | 1.000 (0.000)
Bootstrap CV | 5033.700 (134.233) | 1.000 (0.000)
Model M3: $n=100,p=500$ | Repro samples | 4.710 (0.558) | 0.995 (0.005)
Bootstrap AIC | 5088.030 (456.021) | 1.000 (0.000)
Bootstrap BIC | 2944.325 (245.670) | 1.000 (0.000)
Bootstrap CV | 6458.345 (570.104) | 1.000 (0.000)

Table 1 compares the performances of the model candidate sets produced by the
proposed repro samples approach and the typical residual bootstrap approaches
with different tuning criteria. We report the average cardinality of the model
candidate sets (Cardinality) and the percentage of simulation cases where the
true model $\tau_{0}$ is included in the model candidate sets (Inclusion of
$\tau_{0}$). From Table 1, we see that the proposed repro samples approach
provides much smaller model candidate sets. For model (M1) where
$n=50,p=1000$, with the repro samples approach we manage to cover the true
model 100% of the times with a model candidate set with size of only 2-3 on
average out of the 1000 repro samples. Even in (M2) and (M3), where weak
signals are present and identifying the true model is notoriously challenging,
the proposed procedure produces candidate set that only consists of around 30
and 5 models on average, but covers the true model 98% and 99.5%of the times
respectively. The bootstrap method, on the other hand, yields 100-200 models
out of 1,000 replicates for (M1), around 1700-4000 models out of 10,000
bootstrap runs for (M2), and around 3000-6500 models out of 100,000 bootstrap
samples for (M3), proportionate to a random search. This table clearly
demonstrate the advantage of our proposed method in getting candidate models.
The size of a candidate model set by our repo sample method is small and
manageable, while those by the corresponding Bootstrap methods are all
unreasonably large, making them inapplicable in practice.

Table 2: Simulation Results of the Model Confidence Sets | Medthod |  Cardinality of $\Gamma^{\tau}_{0.95}$ |  Coverage of $\Gamma^{\tau}_{0.95}$
---|---|---|---
Model M1: $n=50,p=1000$ | Repro samples | 2.180 (0.102) | 1.000 (0.000)
Bootstrap AIC | 165.960 (7.201) | 1.000 (0.000)
Bootstrap BIC | 110.440 (5.423) | 1.000 (0.000)
Bootstrap CV | 207.500 (7.890) | 1.000 (0.000)
Model M2: $n=80,p=150$ | Repro samples | 12.050 (0.708) | 0.955 (0.015)
Bootstrap AIC | 4350.850 (134.000) | 0.995 (0.005)
Bootstrap BIC | 2286.015 ( 77.898) | 0.995 (0.005)
Bootstrap CV | 5033.700 (134.233) | 0.995 (0.005)
Model M3: $n=100,p=500$ | Repro samples | 3.810 (0.316) | 0.985 (0.009)
Bootstrap AIC | 3481.215 (497.552) | 0.935 (0.017)
Bootstrap BIC | 1425.125 (265.383) | 0.905 (0.021)
Bootstrap CV | 4989.195 (618.298) | 0.955 (0.015)

Table 2 reports the average cardinality of the confidence sets obtained using
Algorithm 2 and their coverage of the true model $\tau_{0}$ out of the 200
repetitions. From Table 2, we can see that, for Model (M1), the model
confidence set based on the repro samples approach only contains two models on
average, while the sets by the bootstrap methods have sizes between 110-210.
The coverage of the proposed approach is 100%, because the discrete nature of
the model confidence sets could often preclude an exact confidence level of
95%. For example, if a candidate set only contains two models, excluding one
of them could drop the confidence level well below 95%. For Model (M2) and
(M3), the model confidence sets generated by the bootstrap are impractically
large, containing between 1400-5000 models on average. Even with those many
models, the bootstrap confidence set with AIC and BIC slightly under covers
the true model $\tau_{0}$ for (M3). On the contrary, for (M2) and (M3), the
proposed approach manages to give a confidence set of around just 12 and 4
models on average with a good coverage above 95% under the decaying signal
setting, where the minimum coefficient is only 0.6.

The above results demonstrate that the proposed repro samples method
constructs valid and efficient model confidence sets for the true model
$\tau_{0},$ even under the challenging settings of (M1) - (M3). On the
contrary, the bootstrap method exaggerates the uncertainty of model selection
by producing extremely large model confidence sets, rendering results that are
not very useful in practice.

###  4.2 Inference for regression coefficients accounting for model selection
uncertainty

In this subsection, we evaluate the performance of the proposed repro
confidence set for single regression coefficients and compare it with an
existing state-of-art debiased method. We also evaluate the performance of the
joint confidence set for $\bm{\beta}_{0}^{full}.$

Table 3: Simulation Results of Confidence Sets for Regression Parameters |  | Repro Samples | Debiased LASSO
---|---|---|---
Model | $\beta_{0,i}$ | Coverage | Width | Coverage | Width
M1: $n=50,p=1,000$ | All $\beta_{0,i}$ | 1.000(0.000) | 0.003(0.000) | 0.984(0.000) | 0.247(0.000)
$\beta_{0,i}\neq 0$ | 0.952(0.009) | 0.748(0.006) | 0.357(0.020) | 0.247(0.005)
$\beta_{0,i}=0$ | 1.000(0.000) | 0.001(0.000) | 0.986(0.000) | 0 .247(0.000)
M2: $n=80,p=150$ | All $\beta_{0,i}$ | 0.999(0.000) | 0.038(0.000) | 0.982(0.001) | 0.498(0.000)
$\beta_{0,i}\neq 0$ | 0.975(0.005) | 0.539(0.003) | 0.938(0.008) | 0.497(0.002)
$\beta_{0,i}=0$ | 1.000(0.000) | 0.021(0.000) | 0.984(0.001) | 0.498(0.000)
M3: $n=100,p=500$ | All $\beta_{0,i}$ | 1.000(0.000) | 0.007(0.000) | 0.989(0.000) | 0.429(0.000)
$\beta_{0,i}\neq 0$ | 0.959(0.006) | 0.459(0.002) | 0.868(0.010) | 0.428(0.001)
$\beta_{0,i}=0$ | 1.000(0.000) | 0.002(0.000) | 0.991(0.000) | 0.429(0.000)
Table 4: Performances of Repro Sampling vs Debiased LASSO on Inference of the Nonzero Regression Coefficients |  | Repro Samples | Debiased LASSO
---|---|---|---
Model | $\beta_{0,i}$ | Coverage | Width | Coverage | Width
M1: $n=50,p=1,000$ | $\beta_{0,1}=3$ | 0.970(0.012) | 0.714(0.008) | 0.310(0.033) | 0.244(0.008)
$\beta_{0,2}=2$ | 0.960(0.014) | 0.810(0.010) | 0.440(0.035) | 0.249(0.009)
$\beta_{0,3}=1.5$ | 0.925(0.019) | 0.718(0.009) | 0.320(0.033) | 0.250(0.009)
M2: $n=80,p=150$ | $\beta_{0,1}=2$ | 0.990(0.007) | 0.540(0.006) | 0.960(0.014) | 0.498(0.004)
$\beta_{0,2}=1.5$ | 0.965(0.013) | 0.544(0.006) | 0.915(0.020) | 0.498(0.005)
$\beta_{0,3}=1$ | 0.965(0.013) | 0.538(0.007) | 0.935(0.017) | 0.494(0.005)
$\beta_{0,4}=0.8$ | 0.980(0.010) | 0.540(0.006) | 0.930(0.018) | 0.494(0.004)
$\beta_{0,5}=0.6$ | 0.975(0.011) | 0.533(0.007) | 0.950(0.015) | 0.499(0.004)
M3: $n=100,p=500$ | $\beta_{0,1}=3$ | 0.980(0.010) | 0.458(0.005) | 0.870(0.024) | 0.427(0.003)
$\beta_{0,2}=2$ | 0.950(0.015) | 0.463(0.005) | 0.870(0.024) | 0.428(0.003)
$\beta_{0,3}=1.5$ | 0.965(0.013) | 0.461(0.005) | 0.865(0.024) | 0.427(0.004)
$\beta_{0,4}=1$ | 0.960(0.014) | 0.465(0.005) | 0.850(0.025) | 0.431(0.004)
$\beta_{0,5}=0.8$ | 0.945(0.016) | 0.471(0.006) | 0.890(0.022) | 0.431(0.003)
$\beta_{0,6}=0.6$ | 0.955(0.015) | 0.436(0.005) | 0.865(0.024) | 0.424(0.003)

In Table 3, we compare the performance of the repro samples approach in terms
of the constructed confidence sets for single regression coefficients with the
stat-of-the-art debiased/de-sparsify approach [35]. We notice that although
the debiased LASSO could achieve the desired average coverage rate, it
significantly under covers the non-zero coefficients (signals), both for (M2)
with $n=100,p=500$, and more significantly for the extreme high-dimensional
setting (M3) where $n=50,p=1000.$ The proposed repro samples approach on the
other hand successfully achieves the desired coverage for the signals for all
the three models. Moreover, compared to the debiased approach, the confidence
sets for the coefficients whose true values are zero produced by the repro
samples approach appear to be much narrower, making the overall average width
of the proposed confidence sets much smaller. After all, it is evident that
when it comes to construct confidence sets for single coefficients, the repro
samples approach could achieve a better coverage with narrower intervals.

To further investigate the differences in between the performance of the
proposed approach and that of the debiased LASSO, we compare in Table 4 the
coverage rates and widths of the confidence sets produced by both approach for
each of the nonzero regression coefficients. We observe that the repro samples
approach uniformly achieves the desired coverage regardless of the signal
strength, as opposed to the fact that the debiased LASSO uniformly under
covers the truths for all signals in all three models except for two
coefficients in Model (M2). As expected, the under coverage issues of the
debaised LASSO approach is more serious when $p/n$ is larger, since the second
approximation is more difficult. For Model (M1) with $n=50,p=1000$, the
coverage rate of the debiased LASSO is only around 30%-45%, and for Model (M3)
with $n=100,p=500$ the coverage rate is around 85%-89%. When $p/n$ is relative
smaller as in Model (M2) with $n=80,p=150$, the coverage of debiased lasso
improves to over 91% although the under coverage issue still lingers. In terms
of the width of the confidence sets, when $n=100,p=500$, widths of the repro
samples confidence sets are comparable to those from the debiased LASSO, even
though the coverage of the repro samples confidence sets are $10\%$ higher in
general. For Model (M1) with $n=50,p=1000$, the repro samples interval are
wider while maintaining the desired coverages, fairly indicating a larger
uncertainty than (M3) with $n=100,p=500,$ since we have a smaller sample and
much more parameters to deal with. To sum up, the repro samples approach is
able to cover all the signals with desired coverage rate and correctly
quantify the uncertainty of parameter estimation regardless of the dimension
of the design matrix and signal strength.

The proposed repro samples method is more flexible than existing methods.
Beside to get the model confidence set and the confidence set for single
coefficients, it provides a joint inference for $\bm{\beta}_{0}^{full}$. To
evaluate the performance of the joint confidence set for
$\bm{\beta}_{0}^{full}$ in (3.2), we apply (3.2) on the 200 simulated data
sets for models (M1)-(M3), and summarize the results in Table 5. Evidently,
the proposed confidence set can achieve the desired coverage rate, since it
covers the truth $\bm{\beta}_{0}^{full}$ 94%, 94.5% and 95% of the times for
models (M1), (M2) and (M3), respectively. Moreover, the proposed confidence
set, as opposed to those in [56, 21], has a sparse structure in the sense that
the vast majority of dimensions of the joint confidence set corresponding to
the zero regression coefficients are shrunk to [0,0], as illustrated by Table
5. This is because if variable $X_{i}$ is not in any of the models in the
model candidate set $S^{(d)},$ then any value of $\bm{\beta}$ with nonzero
$\beta_{i}$ will be excluded from the confidence set following from the union
in (3.2). Such sparse confidence sets give researchers two advantages in
practice: (1) the size/volume of the confidence set is substantially smaller,
and therefore it is more informative; (2) it offers a new tool for confidently
and efficiently screening variables. Here the proportions of the confidence
set’s dimensions shrunk to $[0,0]$ are above 98.5% for model (M1) and (M3) and
91.6% for model (M2), demonstrating that the variables left after screening is
much smaller than $n-1,$ which is suggested for the sure independence
screening approach [24, 25].

Table 5: Performances of the Joint Confidence Set for $\bm{\beta}_{0}^{full}$ Model | Coverage Rate | Proportions of Dimensions Shrunk to [0,0]
---|---|---
M1: $n=50,p=1,000$ | 0.940 (0.016) | 0.997 (0.000)
M2: $n=80,p=150$ | 0.945 (0.016) | 0.916 (0.002)
M3: $n=100,p=500$ | 0.950 (0.015) | 0.986 (0.001)

##  5 Technical proofs

In this section, we prove our main results Theorems 1-7. We would like to
point out that Theorems 1 and 2 in Section 2.1 are particularly challenging.
In both cases, we have to control the behaviour of the repro samples
$\mathbf{U}^{*},$ not only in relation to the error term $\mathbf{U},$ but
also in relation to
$(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{X}_{\tau_{0}}\beta_{0}$ for any
$|\tau|\leq|\tau_{0}|,$ within the proximity of which $\mathbf{U}^{*}$ could
possibly lead to $\tau$ instead of $\tau_{0}.$ We also would like to note that
there has not been any finite-sample theories like Theorem 1 in the
literature.

In the following, we first prove Theorems 1 and 2, that is, when either the
number of repro simulations $d$ or the sample size $n$ is large enough, our
candidate model sets $S^{(d)}$ will contain the true model $\tau_{0}$ with
high probability.

###  5.1 Proof of Theorem 1

First we define a similarity measure between two vectors
$\mathbf{v}_{1},\mathbf{v}_{2}$ as the square of cosine of the angles between
$\mathbf{v}_{1}$ and $\mathbf{v}_{2},$ i.e.

$\displaystyle\rho(\mathbf{v}_{1},\mathbf{v}_{2})=\frac{\|\mathbf{H}_{\mathbf{v}_{1}}\mathbf{v}_{2}\|^{2}}{\|\mathbf{v}_{2}\|^{2}}=\frac{(\mathbf{v}_{1}^{\top}\mathbf{v}_{2})^{2}}{\|\mathbf{v}_{1}\|^{2}\|\mathbf{v}_{2}\|^{2}}.$

We therefore use $\rho(\mathbf{u}^{*},\mathbf{u}^{rel})$ to measure the
similarity between a single repro sample $\mathbf{u}^{*}$ and the realization
$\mathbf{u}^{rel}$. Apparently, the closer $\rho(\mathbf{u},\mathbf{u}^{rel})$
is to 1, the smaller the angle between $\mathbf{u}$ and $\mathbf{u}^{rel}.$
Hence we use $\rho(\mathbf{u},\mathbf{u}^{rel})$ to measure the similarity
between $\mathbf{u}$ and $\mathbf{u}^{rel}$.

We then present a technical lemma that derives the probability bound of
obtaining the true model $\tau_{0}$ when the repro sample $\mathbf{u}^{*}$
falls within close proximity of $\mathbf{u}^{rel}$ in that
$\rho(\mathbf{u}^{*},\mathbf{u}^{rel})>1-\gamma^{2}_{2}$ for a small
$\gamma_{2}>0.$ We provide the proof of Lemma 2 in the appendix.

######  Lemma 2.

Suppose $n-|\tau_{0}|>4.$ Let $\mathbf{U}^{*}$ be a random repro sample of
$\mathbf{U}$, such that $\mathbf{U}^{*},\mathbf{U}\sim N(0,\mathbf{I}_{n}),$
and

$\hat{\tau}_{\mathbf{U}^{*}}=\mathop{\rm
arg\,min\,}_{\tau}\left\\{\min_{\bm{\beta}_{\tau},\sigma}\|\mathbf{Y}-\mathbf{X}_{\tau}\bm{\beta}_{\tau}-\sigma\mathbf{U}^{*}\|^{2}+\lambda|\tau|\right\\}.$

Then for any
$0<\gamma^{1/4}_{2}<\min\left\\{\frac{C_{\min}}{|\tau_{0}|^{2}},\frac{C_{\min}}{24\\{2+2(|\tau_{0}|+1)\log(p/2)/n\\}},0.35\right\\}$
such that
$C_{\min}>52{\sqrt{\gamma_{2}}}\left(\frac{\log(p/2)}{n}+\gamma_{2}\right)\sigma_{0}^{2}$
and
$\lambda\in\left[n\gamma_{2}^{1/2}\left\\{2+2(|\tau_{0}|+1)\frac{\log(p/2)}{n}\right\\},n\gamma_{2}^{1/4}\frac{C_{\min}}{6}\right],$

$\displaystyle\mathbb{P}\left\\{\hat{\tau}_{\mathbf{U}^{*}}\neq\tau_{0},\rho(\mathbf{U}^{*},\mathbf{U})>1-\gamma_{2}^{2}\right\\}$
$\displaystyle\leq
3\exp\left\\{-\frac{n}{26\sigma_{0}^{2}}\left[\frac{C_{\min}}{\sqrt{\gamma_{2}}}-52\left(\frac{\log(p/2)}{n}+\gamma_{2}\right)\sigma_{0}^{2}\right]\right\\}+3\exp\left(-\frac{n}{4\gamma_{2}^{1/4}\sigma_{0}^{2}}\right)$
(24) $\displaystyle\hskip
113.81102pt+4(64\gamma_{2})^{\frac{n-|\tau_{0}|-1}{6}}p^{|\tau_{0}|}+2\gamma_{2}^{\frac{n-|\tau_{0}|}{2}-1}(\sqrt{n}p)^{|\tau_{0}|}.$

Unlike existing literature in high-dimensional regime, the results in Lemma 2
do not require any conditions on $C_{\min},$ nor do they even depend on any
conditions necessary for achieving consistent regression parameter estimation.
This is because the probability bound on the right hand side of Eq. (2)
depends on $C_{\min}$ only through ${C_{\min}}/{\sqrt{\gamma_{2}}}$. When the
quantity ${C_{\min}}/{\sqrt{\gamma_{2}}}$ becomes larger, the probability
bound becomes smaller. Therefore no matter how small ${C_{\min}}$ is, as long
as ${C_{\min}>0},$ the quantity ${C_{\min}}/{\sqrt{\gamma_{2}}}$ can be
arbitrarily large when $\gamma_{2}$ approaches 0. Consequently, however small
the separation between the true model $\tau_{0}$ and the alternative models
is, we can always recover $\tau_{0}$ with high probability with a repro sample
$\mathbf{U}^{*}$ that is close to $\mathbf{U}.$

By the finite probability bound obtained in the above lemma, when $\gamma_{2}$
goes to 0, that is, $\mathbf{U}^{*}$ proximate $\mathbf{U}$ more closely, the
probability of $\hat{\tau}_{\mathbf{U}^{*}}\neq\tau_{0}$ goes to 0 for any
finite $n$ and $p$. This indicates that we do not need $\mathbf{U}^{*}$ to hit
$\mathbf{U}$ exactly, rather we would only need $\mathbf{U}^{*}$ to be in a
neighbourhood of $\mathbf{U}$ in order to recover $\tau_{0}$ with high
probability. Additionally we observe that as the sample size $n$ goes larger,
the probability bounds in Lemma 2 decay exponentially. Therefore, for larger
sample, the estimation $\hat{\tau}_{\mathbf{U}^{*}}=\tau_{0}$ with large
probability even for a large $\gamma_{2}.$ As a result, the neighbourhood of
$\mathbf{u}^{rel}$, within which $\mathbf{U}$ yields
$\hat{\tau}_{\mathbf{U}^{*}}=\tau_{0}$ with high probability, will expand as
the sample size $n$ grows larger.

As Lemma 2 shows the probability of a single sample $\mathbf{U}$ being close
to $\mathbf{U}^{*}$, in the following Lemma 3, we develop the probability
bound for at least one of the $d$ independent samples of $\mathbf{U}$ being
close to $\mathbf{U}^{*}$. This probability bound then implies a finite-sample
probability bounds of $\tau_{0}$ not included in the candidate set $S^{(d)}$
constructed by Algorithm 1. The proof of Lemma 3 is deferred to the appendix.

######  Lemma 3.

Suppose $n-|\tau_{0}|>4.$ Then for any

$0<\gamma^{1/4}_{2}<\allowbreak\min\bigg{\\{}\frac{C_{\min}}{|\tau_{0}|^{2}},\allowbreak\frac{C_{\min}}{24\\{2+2(|\tau_{0}|+1)\log(p/2)/n\\}},0.35\bigg{\\}},$

such that
$C_{\min}>52{\sqrt{\gamma_{2}}}\left(\frac{\log(p/2)}{n}+\gamma_{2}\right)\sigma_{0}^{2},$
and
$\lambda\in\bigg{[}n\gamma_{2}^{1/2}\left\\{2+2(|\tau_{0}|+1)\frac{\log(p/2)}{n}\right\\},\allowbreak
n\gamma_{2}^{1/4}\frac{C_{\min}}{6}\bigg{]},$ the finite-sample probability
bound that the true model is not covered by the model candidates set
$S^{(d)}$, obtained by Algorithm 1 with the objective function (5), is as
follows,

$\displaystyle\mathbb{P}_{({\cal U}^{d},\mathbf{Y})}(\tau_{0}\notin
S^{(d)})\leq
3\exp\left\\{-\frac{n}{26\sigma_{0}^{2}}\left[\frac{C_{\min}}{\sqrt{\gamma_{2}}}-52\left(\frac{\log(p/2)}{n}+\gamma_{2}\right)\sigma_{0}^{2}\right]\right\\}$
$\displaystyle\quad+3\exp\left(-\frac{n}{4\gamma_{2}^{1/4}\sigma_{0}^{2}}\right)+4(64\gamma_{2})^{\frac{n-|\tau_{0}|-1}{6}}p^{|\tau_{0}|}+2\gamma_{2}^{\frac{n-|\tau_{0}|}{2}-1}(\sqrt{n}p)^{|\tau_{0}|}+\left(1-\frac{\gamma_{2}^{n-1}}{n-1}\right)^{d}.$
(25)

We are now to present the proof of Theorem 1

###### Proof of Theorem 1:.

The first four terms of $\eqref{eq:prob_bound_penalize_finite}$ go to 0 as
$\gamma_{2}$ goes to $0.$ Therefore for any $\delta>0,$ there exists a
$\gamma_{\delta}>0,$ such that when $\gamma_{2}=\gamma_{\delta},$ sum of the
first four terms of $\eqref{eq:prob_bound_penalize_finite}$ is smaller than
$\delta,$ which implies the probability bound in (7).

∎

###  5.2 Proof of Theorem 2

Similar to the last section, we first introduce a key lemma.

######  Lemma 4.

For any finite $n$ and $p$, if

$\frac{\lambda}{n}\in\left[\frac{3\sigma_{0}^{2}(|\tau_{0}|+1)(\log(p-|\tau_{0}|)+\log(|\tau_{0}|)+\frac{2}{3})}{n}+t,\frac{(1-\gamma_{1}^{2})C_{\min}}{6},\right],$

a finite probability bound that the true model is not covered by the model
candidates set $S^{(d)}$, obtained by Algorithm 1 with the objective function
(5), is as follows,

$\displaystyle\mathbb{P}_{({\cal U}^{d},\mathbf{Y})}(\tau_{0}\notin
S^{(d)})\leq
L(\gamma_{1})+3\exp\left(-\frac{nt}{3\sigma_{0}^{2}}\right)+\left[2\\{\arccos(\gamma_{1})\\}^{n-|\tau_{0}|-1}p^{|\tau_{0}|}\right]^{d},$
(26)

where

$\displaystyle
L(\gamma_{1})=6\exp\left[-\frac{n}{18\sigma_{0}^{2}}\left\\{(1-\gamma_{1}^{2})C_{\min}-36\frac{\log
p}{n}\sigma_{0}^{2}\right\\}\right],$

and $\cos(0.3\pi)<\gamma_{1}<1$ is any real number.

Lemma 4 aims to offer insights on asymptotic property of the candidate set
$S^{(d)},$ therefore it gives a different probability bound than Lemma 3. The
interpretation is that for any fixed $d,$ the probability of $\tau_{0}\notin
S^{(d)}$ is $O(e^{-n})$ under the conditions in Theorem 2. This provides us
the insight that for large samples, we actually do not need a extremely large
number of repro samples in order to recover the true model in the candidate
set $S^{(d)}.$

To explain the intuition behind the probability bound in Lemma 4, we denote
the angle between the repro sample $\mathbf{U}^{*}$ and
$(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{X}_{\tau_{0}}\bm{\beta}_{0}$ as
$\gamma_{1}^{\tau}.$ If $\gamma_{1}^{\tau}\geq\gamma_{1}$ for all
$|\tau|\leq|\tau_{0}|,$ then the probability of
$\hat{\tau}_{\mathbf{U}^{*}}\neq\tau_{0}$ is bounded by the first two terms of
(26). The reason that we want to bound $\mathbf{U}^{*}$ away from
$(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{X}_{\tau_{0}}\bm{\beta}_{0}$ is that
when
$\mathbf{U}^{*}\approx(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{X}_{\tau_{0}}\bm{\beta}_{0},$
$\mathbf{X}_{\tau}$ will explain $\mathbf{Y}-\mathbf{U}^{*}$ as well as
$\mathbf{X}_{\tau_{0}},$ possibly leading to
$\hat{\tau}_{\mathbf{U}^{*}}=\tau\neq\tau_{0}.$ The last term of (26) is
derived from the probability bound that $\gamma_{1}^{\tau}\leq\gamma_{1}$ for
some $|\tau|\leq|\tau_{0}|$ for all the $d$ copies of repro samples
$\mathbf{U}^{*}.$ Therefore, all the three terms together give a probability
bound for $\tau_{0}\notin S^{(d)}.$

We now present the proof of Theorem 2.

###### Proof of Theorem 2.

By Lemma 4, we obtain (2) by making $\gamma^{2}_{1}=0.7$. The lower bound for
$\frac{\lambda}{n}$ is simplified by applying the inequality
$\log(|\tau_{0}|)+\log(p-|\tau_{0}|)\leq 2\log(p/2).$ ∎

We then proceed to the proof of Theorem 3-5, showing the validity of the
constructed confidence set.

###  5.3 Proof of Theorems 3\- 5

###### Proof of Theorem 3.

The proof of Theorem 3 is a direct consequence of our repro samples idea.

$\displaystyle\mathbb{P}(\tau_{0}\in\Gamma^{\tau}_{\alpha}(\mathbf{Y}))$
$\displaystyle\geq\mathbb{P}(T(\mathbf{U},\bm{\theta})\in
B_{\alpha}(\tau,\mathbf{W}(\mathbf{U},\bm{\theta}))$
$\displaystyle=\mathbb{E}\\{\mathbb{P}(T(\mathbf{U},\bm{\theta})\in
B_{\alpha}(\tau,\mathbf{W}(\mathbf{U},\bm{\theta}))|\mathbf{W}(\mathbf{U},\bm{\theta}))\\}\geq\alpha$

∎

###### Proof of Theorem 4.

First, for a given $\tau$, the distribution of
$(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}/\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}\|=(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{Y}_{\bm{\theta}}/\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{Y}_{\bm{\theta}}\|$
is free of $(\bm{\beta}_{\tau},\sigma).$ Therefore
$(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{Y}_{\bm{\theta}}/\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{Y}_{\bm{\theta}}\|$
is ancillary for $(\bm{\beta}_{\tau},\sigma).$ Because
$\widetilde{\mathbf{W}}(\mathbf{Y}_{\bm{\theta}},\tau)$ is minimal sufficient
for $(\bm{\beta}_{\tau},\sigma),$ then by Basu’s theorem
$\widetilde{\mathbf{W}}(\mathbf{Y}_{\bm{\theta}},\tau)$ is independent of
$(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}/\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}\|.$
Apparently,
$\mathbf{A}_{\bm{\theta}}(\mathbf{U})=\widetilde{\mathbf{A}}_{\tau}(\mathbf{Y}_{\bm{\theta}})$
and $b_{\bm{\theta}}(\mathbf{U})=\tilde{b}_{\tau}(\mathbf{Y}_{\bm{\theta}})$
are independent. It then follows that
$\widetilde{\mathbf{A}}_{\tau}(\mathbf{Y}_{\bm{\theta}}),$
$\tilde{b}_{\tau}(\mathbf{Y}_{\bm{\theta}})$ and
$(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}/\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}\|$
are mutually independent. As a result, we conclude that the conditional
distribution

$\displaystyle\left\\{\mathbf{Y}_{\bm{\theta}}|\widetilde{\mathbf{W}}(\mathbf{Y}_{\bm{\theta}},\tau)=(\mathbf{a}_{obs},b_{obs})\right\\}\sim\left\\{\mathbf{a}_{obs}+b_{obs}\frac{(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}}{\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}\|}\right\\}\sim\mathbf{Y}^{*},$

where
$\mathbf{Y}^{*}=\left\\{\mathbf{a}_{obs}+b_{obs}\frac{(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}^{*}}{\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}^{*}\|}\right\\}$
and $\mathbf{U}^{*}\sim\mathbf{U},$ is free of $(\bm{\beta}_{\tau},\sigma)$
for any $\mathbf{a}_{obs},b_{obs}.$ Then the conditional probability in (14)
is free of $(\bm{\beta}_{\tau},\sigma),$ hence the Borel set
$B_{\alpha}(\tau,\mathbf{w})$ defined (15) is also free of
$(\bm{\beta}_{\tau},\sigma).$ Moreover, it follows from (15) that

$\displaystyle{\mathbb{P}}_{\mathbf{Y}_{\bm{\theta}}|\mathbf{w}}\left\\{\tilde{T}(\mathbf{Y}_{\theta},\tau)\in
B_{\alpha}(\tau,\mathbf{w})\big{|}\widetilde{\mathbf{W}}(\mathbf{Y}_{\bm{\theta}},\tau)=\mathbf{w}\right\\}$
$\displaystyle=\sum_{\tau^{*}\in
B_{\alpha}(\tau,\mathbf{w})}p_{(\mathbf{w},\tau)}(\tau^{*})$
$\displaystyle=1-\sum_{\tau^{*}\not\in
B_{\alpha}(\tau,\mathbf{w})}p_{(\mathbf{w},\tau)}(\tau^{*})\geq\alpha,$

which proves (16). Then following from (16), (17) and Theorem 3,

$\displaystyle\mathbb{P}_{\mathbf{Y}}\left\\{\tau_{0}\not\in{\bar{\Gamma}}^{\tau}_{\alpha}(\mathbf{Y})\right\\}\leq\mathbb{P}_{\mathbf{Y}}\left\\{\tau_{0}\not\in{\Gamma}^{\tau}_{\alpha}(\mathbf{Y})\right\\}+\mathbb{P}_{\mathbf{Y}}\left\\{\tau_{0}\not\in
S^{(d)}\right\\}\leq 1-\alpha+\mathbb{P}_{\mathbf{Y}}\left\\{\tau_{0}\not\in
S^{(d)}\right\\}.$ (27)

Then Theorem 4 follows from the above and Theorem 1 for some
$c_{1}<-\log\left(1-\frac{\gamma_{\delta}^{n-1}}{n-1}\right).$ ∎

###### proof of Theorem 5.

Under the conditions in Theorem 2, we make the constant $c_{2}$

$\displaystyle
0<c_{2}<\min\left\\{\frac{1}{18\sigma_{0}^{2}}\left(0.3C_{\min}-36\frac{\log(p+1)}{n}\sigma_{0}^{2}\right),\frac{t}{3\sigma_{0}^{2}},d\left(0.23-\frac{|\tau_{0}|\log(p)+2}{n}\right)\right\\}.$
(28)

Then Theorem 5 follows from (27) and Theorem 2. ∎

Next, we present the proofs of Theorems 6 and 7, showing the validity of the
inference for any subset of regression coefficients, both in finite samples
and asymptotically.

###  5.4 Proofs of Theorems 6 and 7

###### Proof of Theorem 6.

We first write

$\displaystyle\mathbb{P}_{\mathbf{Y}}\\{\bm{\beta}_{0,\Lambda}\not\in\Gamma^{\bm{\beta}_{\Lambda}}_{\alpha}(\mathbf{Y})\\}$
$\displaystyle=\mathbb{P}_{\mathbf{Y}}\\{\bm{\beta}_{0,\Lambda}\not\in\Gamma^{\bm{\beta}_{\Lambda}}_{\alpha}(\mathbf{Y}),\tau_{0}\in
S^{(d)}\\}+\mathbb{P}_{\mathbf{Y}}\\{\beta_{0,\Lambda}\not\in\Gamma^{\beta_{\Lambda}}_{\alpha}(\mathbf{Y}),\tau_{0}\not\in
S^{(d)}\\}.$ (29)

Then let $\bm{\eta}_{0,\Lambda}=(\tau_{0},\beta_{0,\Lambda}),$

$\displaystyle
P_{\mathbf{Y}}\\{\bm{\beta}_{0,\Lambda}\not\in\Gamma^{\bm{\beta}_{\Lambda}}_{\alpha}(\mathbf{Y}),\tau_{0}\in
S^{(d)}\\}\leq\mathbb{P}_{\mathbf{U}}\\{\tilde{T}(\mathbf{Y},\bm{\eta}_{0,\Lambda})\in
B_{\bm{\eta}_{0,\Lambda}}(\alpha)\\}=1-\alpha.$ (30)

Therefore, from the above and Theorem 1,

$\displaystyle\mathbb{P}_{\mathbf{Y}}\\{\beta_{0,\Lambda}\not\in\Gamma^{\beta_{\Lambda}}_{\alpha}(\mathbf{Y}),\tau_{0}\not\in
S^{(d)}\\}\leq\mathbb{P}_{\mathbf{Y}}\\{\tau_{0}\not\in
S^{(d)}\\}=\delta+o(e^{-c_{1}d}),$

for some $c_{1}<-\log\left(1-\frac{\gamma_{\delta}^{n-1}}{n-1}\right).$
Theorem 6 then follows immediately. ∎

###### Proof of Theorem 7.

Define the constant $c_{2}$ as in (28), then it follows from Theorem 2 that

$\displaystyle\mathbb{P}_{\mathbf{Y}}\\{\beta_{0,\Lambda}\not\in\Gamma^{\beta_{\Lambda}}_{\alpha}(\mathbf{Y}),\tau_{0}\not\in
S^{(d)}\\}\leq\mathbb{P}_{\mathbf{Y}}\\{\tau_{0}\not\in
S^{(d)}\\}=o(e^{-c_{2}n}).$

Hence Theorem 7 follows from (29), (30) and the above. ∎

##  6 Discussion

In this article, we develop a repro samples method to solve the inferential
problems concerning high-dimensional linear regression models. We propose
novel statistical methods and algorithms, with supporting theories, to
construct performance-guaranteed confidence sets for the unknown underlying
true model, as well as for a single or a subset of regression coefficients.

Our work provides a solution to a challenging open problem in the high-
dimensional statistics literature: how to infer the true model and quantify
model selection uncertainty. Besides, the proposed confidence sets
theoretically guarantee the coverage rate in both finite and large samples, as
opposed to existing literature that usually rely on asymptotic results
requiring the sample size go to infinity, In addition, numerical studies
demonstrate that the proposed methods are valid and efficient, and they
outperform the existing methods.

Furthermore, this paper contributes three key technical innovations.

  1. 1.

We develop a data-driven approach to obtain a efficient model candidate set,
which is able to cover the true model with high probability by including just
a reasonable number of model candidates. Using this model candidate set
effectively addresses the computational issue since it avoids searching the
entire model space. The approach is based on the matching attempt of repro
samples with the observed data, leading to the many-to-one mapping function in
(5). Specifically, this many-to-one mapping tells us that there always exists
a neighbourhood of $\mathbf{U},$ within which a repro copy $\mathbf{U}^{*}$
can help recover the true model with a high probability. With this insight, we
propose a formal procedure and provide supporting theories and numerical
evidence, both of which also help to outline trade-offs among sample size, the
signal strength, and the performance of the model candidate set.

  2. 2.

When making inference for the unknown model, we develop a conditional repro
samples approach to remove the impact of the nuisance parameters
$(\bm{\beta}_{\tau},\sigma^{2})$. This conditional approach works in general
for inference problems beyond the scope of this paper. In particular, let
$\bm{\theta}_{0}=(\bm{\nu}_{0},\bm{\xi}_{0})$, where $\bm{\nu}_{0}$ and
$\bm{\xi}_{0}$ are interested and nuisance parameters, respectively. If we
have a nuclear mapping $T(\mathbf{U},\bm{\theta})$ and a quantity
$\mathbf{W}(\mathbf{U},\bm{\theta}),$ such that the conditional distribution
of $T(\mathbf{U},\bm{\theta})$ given $W(\mathbf{U},\bm{\theta})=\mathbf{w}$ is
free of $\bm{\xi},$ then there exists a Borel set
$B_{\alpha}(\bm{\nu},\mathbf{w})$ free of the nuisance $\bm{\xi}$ such that
$\mathbb{P}(T(\mathbf{U},\bm{\theta})\in
B_{\alpha}(\bm{\nu},w)|W(\mathbf{U},\bm{\theta})=w)\geq\alpha.$ Then similar
to (2.2), we can get a valid level-$\alpha$ confidence set for $\bm{\nu}_{0}$.

  3. 3.

We propose confidence sets both for a single and for any subset of regression
coefficients. In contrast, existing literature only focus on one aspect of
these inference problems. This is because unlike existing approaches, we take
a union of intervals or multi-dimensional ellipsoids based on each model in
the model candidate set. Thereby, our approach takes into account the
uncertainty in model estimation. Not only does it provides the desired
coverage, it also produces confidence sets that are sparse and generally
smaller than the existing methods, including the debiased approach.

Finally, there are several potential directions for extensions of the work.
First, in this paper we considered the Gaussian model (2) with $\mathbf{U}\sim
N(0,1)$. It is possible to extend it to other given distributions or a
generalized linear model setting. In these extensions, the conditional
approach to handle nuisance regression parameters may not directly apply, but
we may use asymptotic arguments to justify their validity and coverage
statements. Moreover, another research direction is how to handle weak
signals. Even though we do not need any condition on the signal strength, the
proposed approaches may demand large computational cost when weak signals are
present. Therefore a natural question is, under limited computational
resource, how to adjust the proposed approaches for weak signals.
Additionally, the identifiability condition or $C_{\min}>0$ ensures that there
is no perfect co-linearity between the true model and an alternative model of
equal size. When there is, then multiple equivalent “true” models may exist.
The proposed approach is still valid to cover one of these “true” models.
Further research, however, is still needed if we would like to have a
confident set to cover all of these equivalent “true” models. These topics are
outside the scope of the current paper and will be reported in other articles.

##  Appendix A An alternative formulation in Section 2

In Section 2, we obtain the candidate sets by solving the objective function
(5). In addition to (5), there is also an almost equivalent form that imposes
a constraint on $|\tau|=\|\bm{\beta}_{\tau}\|_{0}$ other than adding a
regularization term, i.e.

$\displaystyle\min\|\mathbf{y}_{obs}-X_{\tau}\bm{\beta}_{\tau}-\sigma\mathbf{U}^{*}_{b}\|^{2}_{2},\quad\hbox{s.t.
}|\tau|\leq k,$ (31)

where $k$ is a constraint on the model size, playing a similar role as the
$\lambda$ in (5). One can opt to use (31) in the Step 2 of Algorithm 1.
Similarly, We can obtain the following results, where Theorem 8, Lemma 5,
Lemma 6 and Theorem 9 are counterparts of Theorem 1, Lemma 2, Lemma 3 and
Theorem 2 respectively.

######  Theorem 8.

For any $\delta>0$, there exists a constant $\gamma_{\delta}$ such that under
the constraint $|\tau|\leq|\tau_{0}|,$ the finite probability bound that the
true model is not covered by the model candidates set $S^{(d)}$, obtained by
Algorithm 1 with the objective function (31), is as follows,

$\displaystyle\mathbb{P}_{({\cal U}^{d},\mathbf{Y})}(\tau_{0}\notin
S^{(d)})\leq\left\\{1-\frac{(\gamma_{\delta})^{n-1}}{n-1}\right\\}^{d}+\delta.$

Therefore as $d\rightarrow\infty,$ $\mathbb{P}_{({\cal
U}^{d},\mathbf{Y})}(\tau_{0}\notin S^{(d)})\rightarrow 0.$

######  Lemma 5.

Suppose $n-|\tau_{0}|>4.$ Under the constraint $|\tau|\leq|\tau_{0}|,$ let

$\hat{\tau}_{\mathbf{U}}=\mathop{\rm
arg\,min\,}_{\\{\tau||\tau|\leq|\tau_{0}|\\}}\left\\{\min_{\bm{\beta}_{\tau},\sigma}\|\mathbf{Y}-\mathbf{X}_{\tau}\bm{\beta}_{\tau}-\sigma\mathbf{U}\|^{2}\right\\}.$

Then for any $0<\gamma_{2}<1/64$ such that
$C_{\min}>24\sqrt{\gamma}_{2}\left(\frac{\log(p/2)}{n}+\gamma_{2}\right)\sigma_{0}^{2},$

$\displaystyle\mathbb{P}\left\\{\hat{\tau}_{\mathbf{U}}\neq\tau_{0},\rho(\mathbf{U}^{*},\mathbf{U})>1-\gamma_{2}^{2}\right\\}$
$\displaystyle\leq
3\exp\left\\{-\frac{n}{12\sigma_{0}^{2}}\left[\frac{C_{\min}}{\sqrt{\gamma}_{2}}-24\left(\frac{\log(p/2)}{n}+\gamma_{2}\right)\sigma_{0}^{2}\right]\right\\}+4(64\gamma_{2})^{\frac{n-|\tau_{0}|-1}{6}}p^{|\tau_{0}|}$
$\displaystyle\hskip
113.81102pt+2\gamma_{2}^{\frac{n-|\tau_{0}|}{2}-1}(\sqrt{n}p)^{|\tau_{0}|}.$

######  Lemma 6.

Suppose $n-|\tau_{0}|>4.$ Then for any $0<\gamma_{2}<1/64$ such that
$C_{\min}>24\sqrt{\gamma}_{2}\left(\frac{\log(p/2)}{n}+\gamma_{2}\right)\sigma_{0}^{2},$
the finite probability bound that the true model is not covered by the model
candidates set $S^{(d)}$, obtained by Algorithm 1 with the objective function
(31), is as follows,

$\displaystyle\mathbb{P}_{({\cal U}^{d},\mathbf{Y})}(\tau_{0}\notin S^{(d)})$
$\displaystyle\leq
3\exp\left\\{-\frac{n}{12\sigma_{0}^{2}}\left[\frac{C_{\min}}{\sqrt{\gamma}_{2}}-24\left(\frac{\log(p/2)}{n}+\gamma_{2}\right)\sigma_{0}^{2}\right]\right\\}$
$\displaystyle+4(64\gamma_{2})^{\frac{n-|\tau_{0}|-1}{6}}p^{|\tau_{0}|}+2\gamma_{2}^{\frac{n-|\tau_{0}|}{2}-1}(\sqrt{n}p)^{|\tau_{0}|}+\left(1-\frac{\gamma_{2}^{n-1}}{n-1}\right)^{d}.$
(32)

######  Theorem 9.

Under the constraint $|\tau|\leq|\tau_{0}|$, the probability bound that the
true model is not covered by the model candidates set $S^{(d)}$, obtained by
Algorithm 1 with the objective function (31) for any finite $d$ is as follows,

$\displaystyle\mathbb{P}_{({\cal U}^{d},\mathbf{Y})}(\tau_{0}\notin
S^{(d)})\leq
6\exp\left[-\frac{n}{18\sigma_{0}^{2}}\\{0.3C_{\min}-36\frac{\log(p)}{n}\sigma_{0}^{2}\\}\right]$
$\displaystyle+\exp\left\\{-nd\left(0.23-\frac{|\tau_{0}|\log(p)+2}{n}\right)\right\\}.$
(33)

Therefore $\mathbb{P}_{({\cal U}^{d},\mathbf{Y})}(\tau_{0}\notin
S^{(d)})\rightarrow 0$ for any $d$ as $n\rightarrow\infty$ , if
$\frac{|\tau_{0}|\log(p)}{n}<0.23$ and
$C_{\min}>120\frac{\log(p+1)}{n}\sigma_{0}^{2}$ when $n$ is large enough.

##  Appendix B Technical Proofs

###  B.1 Proof of Lemma 1

###### Proof of Lemma 1.

By the definition (3), there exist a $\beta_{0}$ and a $\sigma_{0}$ such that
$X_{\tau_{0}}\beta_{0}={\bf y}_{obs}-\sigma_{0}\mathbf{u}^{rel}$. Since

$\displaystyle 0$ $\displaystyle\leq$
$\displaystyle\underset{\tau,\bm{\beta}_{\tau},\sigma}{\min}\|\mathbf{y}_{obs}-\mathbf{X}_{\tau}\bm{\beta}_{\tau}-\sigma\mathbf{u}^{rel}\|^{2}_{2}\leq\|\mathbf{y}_{obs}-X_{\tau_{0}}\beta_{0}-\sigma_{0}\mathbf{u}^{rel}\|^{2}=0,$

it follows that

$\displaystyle\underset{\tau,\bm{\beta}_{\tau},\sigma}{\min}\|\mathbf{y}_{obs}-\mathbf{X}_{\tau}\bm{\beta}_{\tau}-\sigma\mathbf{u}^{rel}\|^{2}_{2}=0.$

Now, let

$(\tilde{\tau},\beta_{\tilde{\tau}},\sigma_{\tilde{\tau}})=\underset{\tau,\bm{\beta}_{\tau},\sigma}{\rm
argmin}\big{\\{}\lambda|\tau|+\|\mathbf{y}_{obs}-\mathbf{X}_{\tau}\bm{\beta}_{\tau}-\sigma\mathbf{u}^{rel}\|^{2}_{2}\big{\\}}.$
(34)

We show below that
$\|\mathbf{y}_{obs}-X_{\tilde{\tau}}\beta_{\tilde{\tau}}-\sigma_{\tilde{\tau}}\mathbf{u}^{rel}\|^{2}_{2}=0$
using the “proof by contradiction” method.

First, we show that, if
$\|\mathbf{y}_{obs}-X_{\tilde{\tau}}\beta_{\tilde{\tau}}-\sigma_{\tilde{\tau}}\mathbf{u}^{rel}\|^{2}_{2}\neq
0,$ then size of $\tilde{\tau}$ must be smaller than $\tau_{0}$, i.e.
$|\tilde{\tau}|<|\tau_{0}|.$ This is because otherwise if
$|\tilde{\tau}|\geq|\tau_{0}|,$ then
$\lambda|\tilde{\tau}|+\|\mathbf{y}_{obs}-X_{\tilde{\tau}}\beta_{\tilde{\tau}}-\sigma_{\tilde{\tau}}\mathbf{u}^{rel}\|^{2}_{2}>\lambda|\tilde{\tau}|\geq\lambda|\tau_{0}|=\lambda|\tau_{0}|+\|\mathbf{y}_{obs}-X_{\tau_{0}}\beta_{\tau_{0}}-\sigma_{0}\mathbf{u}^{rel}\|^{2}_{2},$
which contradicts with (34).

Now, with the triplet
$(\tilde{\tau},\beta_{\tilde{\tau}},\sigma_{\tilde{\tau}})$ defined in (34)
and $|\tilde{\tau}|<|\tau_{0}|<n$, we have for the given $\tilde{\tau}$,

$\displaystyle\|\mathbf{y}_{obs}-\mathbf{X}_{\tilde{\tau}}\bm{\beta}_{\tilde{\tau}}-\sigma_{\tilde{\tau}}\mathbf{u}^{rel}\|\geq\|(\mathbf{I}-\mathbf{H}_{\tilde{\tau},\mathbf{u}^{rel}})\mathbf{y}_{obs}\|,$

where
$\mathbf{H}_{\tilde{\tau},\mathbf{u}^{rel}}=\mathbf{X}_{\tilde{\tau},\mathbf{u}^{rel}}(\mathbf{X}_{\tilde{\tau},\mathbf{u}^{rel}}^{\top}\mathbf{X}_{\tilde{\tau},\mathbf{u}^{rel}})^{-1}\mathbf{X}_{\tilde{\tau},\mathbf{u}^{rel}}^{\top}$
with
$\mathbf{X}_{\tilde{\tau},\mathbf{u}^{rel}}=(\mathbf{X}_{\tilde{\tau}},\mathbf{u}^{rel})$
is the projection matrix to the space expanded by $\mathbf{X}_{\tilde{\tau}}$
and $\mathbf{u}^{rel}$. It follows that

$\displaystyle\|\mathbf{y}_{obs}-\mathbf{X}_{\tilde{\tau}}\bm{\beta}_{\tilde{\tau}}-\sigma_{\tilde{\tau}}\mathbf{u}^{rel}\|\geq\|(\mathbf{I}-\mathbf{H}_{\tilde{\tau},\mathbf{u}^{rel}})\mathbf{y}_{obs}\|=\|(\mathbf{I}-\mathbf{H}_{\tilde{\tau},\mathbf{u}^{rel}})\mathbf{X}_{0}\bm{\beta}_{0}\|,$
(35)

where the equality holds because $\mathbf{u}^{rel}$ is orthogonal to
$(\mathbf{I}-\mathbf{H}_{\tilde{\tau},\mathbf{u}^{rel}})$.

By (35), the definitions of $\gamma^{2}_{(\mathbf{u}^{rel},\tau_{0})}$ and
$C_{\min}$, and under the condition that $0<\lambda\leq
n\\{1-\gamma^{2}_{(\mathbf{u}^{rel},\tau_{0})}\\}C_{\min},$

$\displaystyle\|\mathbf{y}_{obs}-\mathbf{X}_{\tilde{\tau}}\bm{\beta}_{\tilde{\tau}}$
$\displaystyle-\sigma_{\tilde{\tau}}\mathbf{u}^{rel}\|^{2}+\lambda|\tilde{\tau}|\geq\|(\mathbf{I}-\mathbf{H}_{\tilde{\tau},\mathbf{u}^{rel}})\mathbf{X}_{0}\bm{\beta}_{0}\|^{2}+\lambda|\tilde{\tau}|$
$\displaystyle\geq\left\\{1-\gamma^{2}_{(\mathbf{u}^{rel},\tau_{0})}\right\\}\|(\mathbf{I}-\mathbf{H}_{\tilde{\tau}})\mathbf{X}_{0}\bm{\beta}_{0}\|^{2}+\lambda|\tilde{\tau}|$
$\displaystyle\geq\left\\{1-\gamma^{2}_{(\mathbf{u}^{rel},\tau_{0})}\right\\}\,n|\tau_{0}\setminus\tilde{\tau}|C_{\min}+\lambda|\tau_{0}|-\lambda|\tau_{0}\setminus\tilde{\tau}|$
$\displaystyle\geq\lambda|\tau_{0}|=\lambda|\tau_{0}|+\|\mathbf{y}_{obs}-\mathbf{X}_{\tau_{0}}\bm{\beta}_{0}-\sigma_{0}\mathbf{u}^{rel}\|^{2}_{2},$

which contradicts with (34). Thus,
$\|\mathbf{y}_{obs}-\mathbf{X}_{\tilde{\tau}}\bm{\beta}_{\tilde{\tau}}-\sigma_{\tilde{\tau}}\mathbf{u}^{rel}\|^{2}_{2}\neq
0$ does not hold and we only have
$\|\mathbf{y}_{obs}-\mathbf{X}_{\tilde{\tau}}\bm{\beta}_{\tilde{\tau}}-\sigma_{\tilde{\tau}}\mathbf{u}^{rel}\|^{2}_{2}=0.$
Because
$\mathbf{u}^{rel}\not\in\operatorname{span}(\mathbf{X}_{0},\mathbf{X}_{\tau})$
for any $\tau$ with $|\tau|\leq|\tau_{0}|,$ by definition (3), we have
$\tilde{\tau}=\tau_{0}$ and thus the conclusion of the lemma follows.

∎

###  B.2 Proof of Lemma 2 and Lemma 5

Before we proceed to the proof of Lemma 2 and Lemma 5, we first provide a few
technical lemmas that can facilitate the proof.

######  Lemma 7.

For any $\tau$ and $\mathbf{u}^{*}$,

$\displaystyle\mathbf{I}-\mathbf{H}_{\tau,\mathbf{u}^{*}}=I-\mathbf{H}_{\tau}-\mathbf{O}_{\tau^{\bot}\mathbf{u}^{*}},$

where
$H_{\tau,\mathbf{U}^{*}}=\begin{pmatrix}X_{\tau}&\mathbf{u}^{*}\end{pmatrix}\begin{pmatrix}X_{\tau}^{\top}X_{\tau}&X_{\tau}^{\top}\mathbf{u}^{*}\\\
(\mathbf{u}^{*})^{\top}X_{\tau}&(\mathbf{u}^{*})^{\top}\mathbf{u}^{*}\end{pmatrix}^{-1}\begin{pmatrix}X_{\tau}^{\top}\\\
(\mathbf{u}^{*})^{\top}\end{pmatrix}$ is the projection matrix on the space
spanned by $(X_{\tau},\mathbf{u}^{*})$ and
$\mathbf{O}_{\tau^{\bot}\mathbf{u}^{*}}=\frac{(I-\mathbf{H}_{\tau})\mathbf{u}^{*}(\mathbf{u}^{*})^{\top}(I-\mathbf{H}_{\tau})}{(\mathbf{u}^{*})^{\top}(I-\mathbf{H}_{\tau})\mathbf{u}^{*}}$
is the projection matrix on the space spanned by
$(I-\mathbf{H}_{\tau})\mathbf{u}^{*}.$

###### Proof.

By a direct calculation, we have

$\displaystyle\mathbf{I}-\mathbf{H}_{\tau,\mathbf{U}^{*}}=I-\begin{pmatrix}X_{\tau}&\mathbf{u}^{*}\end{pmatrix}\begin{pmatrix}X_{\tau}^{\top}X_{\tau}&X_{\tau}^{\top}\mathbf{u}^{*}\\\
(\mathbf{u}^{*})^{\top}X_{\tau}&(\mathbf{u}^{*})^{\top}\mathbf{u}^{*}\end{pmatrix}^{-1}\begin{pmatrix}X_{\tau}^{\top}\\\
(\mathbf{u}^{*})^{\top}\end{pmatrix}$
$\displaystyle=I-\begin{pmatrix}X_{\tau}&\mathbf{u}^{*}\end{pmatrix}\begin{pmatrix}(X_{\tau}^{\top}X_{\tau})^{-1}+\frac{(X_{\tau}^{\top}X_{\tau})^{-1}X_{\tau}^{\top}\mathbf{u}^{*}(\mathbf{u}^{*})^{\top}X_{\tau}(X_{\tau}^{\top}X_{\tau})^{-1}}{(\mathbf{u}^{*})^{\top}(I-\mathbf{H}_{\tau})\mathbf{u}^{*}}&-\frac{(X_{\tau}^{\top}X_{\tau})^{-1}X_{\tau}^{\top}\mathbf{u}^{*}}{(\mathbf{u}^{*})^{\top}(I-\mathbf{H}_{\tau})\mathbf{u}^{*}}\\\
-\frac{(\mathbf{u}^{*})^{\top}X_{\tau}(X_{\tau}^{\top}X_{\tau})^{-1}}{(\mathbf{u}^{*})^{\top}(I-\mathbf{H}_{\tau})\mathbf{u}^{*}}&\frac{1}{(\mathbf{u}^{*})^{\top}(I-\mathbf{H}_{\tau})\mathbf{u}^{*}}\end{pmatrix}\begin{pmatrix}X_{\tau}^{\top}\\\
(\mathbf{u}^{*})^{\top}\end{pmatrix}$
$\displaystyle=I-\begin{pmatrix}X_{\tau}&\mathbf{u}^{*}\end{pmatrix}\begin{pmatrix}(X_{\tau}^{\top}X_{\tau})^{-1}X_{\tau}^{\top}+\frac{(X_{\tau}^{\top}X_{\tau})^{-1}X_{\tau}^{\top}\mathbf{u}^{*}(\mathbf{u}^{*})^{\top}X_{\tau}(X_{\tau}^{\top}X_{\tau})^{-1}X_{\tau}^{\top}}{(\mathbf{u}^{*})^{\top}(I-\mathbf{H}_{\tau})\mathbf{u}^{*}}-\frac{(X_{\tau}^{\top}X_{\tau})^{-1}X_{\tau}^{\top}\mathbf{u}^{*}(\mathbf{u}^{*})^{\top}}{(\mathbf{u}^{*})^{\top}(I-\mathbf{H}_{\tau})\mathbf{u}^{*}}\\\
-\frac{(\mathbf{u}^{*})^{\top}X_{\tau}(X_{\tau}^{\top}X_{\tau})^{-1}X_{\tau}^{\top}}{(\mathbf{u}^{*})^{\top}(I-\mathbf{H}_{\tau})\mathbf{u}^{*}}+\frac{(\mathbf{u}^{*})^{\top}}{(\mathbf{u}^{*})^{\top}(I-\mathbf{H}_{\tau})\mathbf{u}^{*}}\end{pmatrix}$
$\displaystyle=I-X_{\tau}(X_{\tau}^{\top}X_{\tau})^{-1}X_{\tau}^{\top}-\frac{X_{\tau}(X_{\tau}^{\top}X_{\tau})^{-1}X_{\tau}^{\top}\mathbf{u}^{*}(\mathbf{u}^{*})^{\top}X_{\tau}(X_{\tau}^{\top}X_{\tau})^{-1}X_{\tau}^{\top}}{(\mathbf{u}^{*})^{\top}(I-\mathbf{H}_{\tau})\mathbf{u}^{*}}$
$\displaystyle+\frac{X_{\tau}(X_{\tau}^{\top}X_{\tau})^{-1}X_{\tau}^{\top}\mathbf{u}^{*}(\mathbf{u}^{*})^{\top}}{(\mathbf{u}^{*})^{\top}(I-\mathbf{H}_{\tau})\mathbf{u}^{*}}+\frac{\mathbf{u}^{*}(\mathbf{u}^{*})^{\top}X_{\tau}(X_{\tau}^{\top}X_{\tau})^{-1}X_{\tau}^{\top}}{(\mathbf{u}^{*})^{\top}(I-\mathbf{H}_{\tau})\mathbf{u}^{*}}-\frac{\mathbf{u}^{*}(\mathbf{u}^{*})^{\top}}{(\mathbf{u}^{*})^{\top}(I-\mathbf{H}_{\tau})\mathbf{u}^{*}}$
$\displaystyle=I-\mathbf{H}_{\tau}-\frac{\mathbf{H}_{\tau}\mathbf{u}^{*}(\mathbf{u}^{*})^{\top}\mathbf{H}_{\tau}}{(\mathbf{u}^{*})^{\top}(I-\mathbf{H}_{\tau})\mathbf{u}^{*}}+\frac{\mathbf{H}_{\tau}\mathbf{u}^{*}(\mathbf{u}^{*})^{\top}}{(\mathbf{u}^{*})^{\top}(I-\mathbf{H}_{\tau})\mathbf{u}^{*}}+\frac{\mathbf{u}^{*}(\mathbf{u}^{*})^{\top}\mathbf{H}_{\tau}}{(\mathbf{u}^{*})^{\top}(I-\mathbf{H}_{\tau})\mathbf{u}^{*}}-\frac{\mathbf{u}^{*}(\mathbf{u}^{*})^{\top}}{(\mathbf{u}^{*})^{\top}(I-\mathbf{H}_{\tau})\mathbf{u}^{*}}$
$\displaystyle=I-\mathbf{H}_{\tau}-\frac{(I-\mathbf{H}_{\tau})\mathbf{u}^{*}(\mathbf{u}^{*})^{\top}(I-\mathbf{H}_{\tau})}{(\mathbf{u}^{*})^{\top}(I-\mathbf{H}_{\tau})\mathbf{u}^{*}}=I-\mathbf{H}_{\tau}-\mathbf{O}_{\tau^{\bot}\mathbf{u}^{*}}.$

∎

Let
$\rho(\mathbf{v}_{1},\mathbf{v}_{2})=\cos^{2}(\mathbf{v}_{1},\mathbf{v}_{2})=\frac{\|\mathbf{H}_{\mathbf{v}_{2}}\mathbf{v}_{1}\|^{2}}{\|\mathbf{v}_{1}\|^{2}}$
be the square of the cosine of the angle between any two $n\times 1$ vectors
$\mathbf{v}_{1}$ and $\mathbf{v}_{2}.$ Further, for any given $\tau$, let
$\rho_{\tau^{\bot}}(\mathbf{v}_{1},\mathbf{v}_{2})=\rho\\{(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{v}_{1},(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{v}_{2}\\}$
be the consine of the angle between
$(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{v}_{1}$ and
$(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{v}_{2}.$

######  Lemma 8.

Suppose $|\tau|<n$. For any $-1\leq\gamma_{1},\gamma_{2}\leq 1,$ if
$\mathbf{U}^{*}\sim\mathbf{U}\sim N(0,\mathbf{I}),$

$\displaystyle\mathbb{P}_{\mathbf{U}^{*}}\left\\{\rho_{\tau^{\bot}}(\mathbf{U}^{*},\mathbf{X}_{0}\bm{\beta}_{0})<\gamma_{1}^{2}\right\\}=\mathbb{P}_{\mathbf{U}}\left\\{\rho_{\tau^{\bot}}(\mathbf{U},\mathbf{X}_{0}\bm{\beta}_{0})<\gamma_{1}^{2}\right\\}>1-2\\{\arccos(\gamma_{1})\\}^{n-|\tau|-1},$

and

$\displaystyle\mathbb{P}_{(\mathbf{U}^{*},\mathbf{U})}\\{\rho(\mathbf{U}^{*},\mathbf{U})>1-\gamma_{2}^{2}\\}>\frac{\gamma_{2}^{n-2}\arcsin(\gamma_{2})}{n-1}.$

Further more $\rho_{\tau^{\bot}}(\mathbf{U}^{*},\mathbf{X}_{0}\bm{\beta}_{0})$
and $\rho(\mathbf{U}^{*},\mathbf{U})$ are independent.

###### Proof.

Let $(\mathbf{I}-\mathbf{H}_{\tau})=\sum_{i=1}^{n-|\tau|}D_{i}D_{i}^{\top}$ be
the eigen decomposition of $(\mathbf{I}-\mathbf{H}_{\tau})$. Denote by
$Z_{i}=D_{i}^{\top}\mathbf{U}$ and
$w_{i}=D_{i}^{\top}\mathbf{X}_{0}\bm{\beta}_{0}$, for $i=1,\ldots,n-|\tau|$.
It follows that $Z_{1},\dots,Z_{n-|\tau|}$ are i.i.d $N(0,1)$ and

$\displaystyle\mathbb{P}_{\mathbf{U}^{*}}\left\\{\rho_{\tau^{\bot}}(\mathbf{U}^{*},\mathbf{X}_{0}\bm{\beta}_{0})<\gamma_{1}^{2}\right\\}=\mathbb{P}_{\mathbf{U}}\left\\{\rho_{\tau^{\bot}}(\mathbf{U},\mathbf{X}_{0}\bm{\beta}_{0})<\gamma_{1}^{2}\right\\}$
$\displaystyle=\mathbb{P}_{\mathbf{U}}\left\\{\frac{\sum_{i=1}^{n-|\tau|}w_{i}Z_{i}}{\sqrt{\sum_{i=1}^{n-|\tau|}w_{i}^{2}}\sqrt{\sum_{i=1}^{n-|\tau|}Z_{i}^{2}}}<\gamma_{1}\right\\}$
$\displaystyle=\mathbb{P}_{\mathbf{U}}\left\\{|\cos(\varphi)|<\gamma_{1}\right\\},$

where $\varphi=\varphi(\mathbf{U})$ (or $\pi-\varphi$) is the angle between
$(Z_{1},\dots,Z_{n-|\tau|})$ and $(w_{1},\dots,w_{n-|\tau|})$ for
$0\leq\varphi\leq\pi$.

We transform the co-ordinates of $Z_{1},\dots,Z_{n-|\tau|}$ into sphere co-
ordinates, with $\varphi$ as the first angle coordinate. It follows from the
Jacobian of the spherical transformation the density function of $\varphi$ is

$f(\varphi)=\sin^{n-|\tau|-2}(\varphi)/c,\quad 0\leq\varphi\leq\pi,$

where
$c=\int_{0}^{\pi}\sin^{n-|\tau|-2}(\varphi)d\varphi=2\int_{0}^{\frac{\pi}{2}}\sin^{n-|\tau|-2}(\varphi)d\varphi=$
is the normalizing constant.

Note that, for $0<\varphi<\pi/2$, we have

$\frac{2}{\pi}\varphi<\sin(\varphi)<\min\\{\varphi,1\\}=\varphi
1_{(0<\varphi<1)}+1_{(1\leq\varphi<\pi/2)},$

where $1_{(\cdot)}$ is an indicator function. It follows that

$\frac{\pi}{2(n-|\tau|-1)}<c<\frac{1}{n-|\tau|-1}+(\frac{\pi}{2}-1)<2.$

Therefore, we have

$\displaystyle\mathbb{P}_{\mathbf{U}}\left\\{|\cos\varphi|<\gamma_{1}\right\\}=\frac{2}{c}\int_{\arccos(\gamma_{1})}^{\pi/2}\sin^{n-|\tau|-2}(s)ds=1-\frac{2}{c}\int_{0}^{{\arccos(\gamma_{1})}}\sin^{n-|\tau|-2}(s)ds$
$\displaystyle\quad>1-\frac{2(n-|\tau|-1)\int_{0}^{{\arccos(\gamma_{1})}}s^{n-|\tau|-2}ds}{\pi}=1-2\\{\arccos(\gamma_{1})\\}^{n-|\tau|-1}.$

Next conditioning on $\mathbf{U}^{*}=\mathbf{u}^{*}$, with similar procedure
as above but replacing $n-|\tau|$ with $n$, we can show that

$\displaystyle\mathbb{P}_{\mathbf{U}}\left\\{{\|(\mathbf{u}^{*})\top\mathbf{U}\|}\big{/}{(\|\mathbf{u}^{*}\|\|\mathbf{U}\|)}>\sqrt{1-\gamma_{2}^{2}}\bigg{|}\mathbf{u}^{*}\right\\}=\mathbb{P}_{\mathbf{U}}\left\\{|\cos(\psi)|>\sqrt{1-\gamma_{2}^{2}}\bigg{|}\mathbf{u}^{*}\right\\}$
$\displaystyle\qquad=\frac{2}{c_{1}}\int_{0}^{\arcsin\gamma_{2}}\sin^{n-2}(s)ds$
(36)
$\displaystyle\qquad>\frac{2}{c_{1}}\int_{0}^{\arcsin\gamma_{2}}(\frac{s\gamma_{2}}{\arcsin\gamma_{2}})^{n-2}ds>\frac{\gamma_{2}^{n-2}\arcsin\gamma_{2}}{n-1},$

where $\psi=\psi(\mathbf{u},^{*}\mathbf{u})$ (or $\pi-\psi$) is the angle
between $\mathbf{u}$ and $\mathbf{u}^{*}$ and the normalizing constant
$c_{1}=\int_{0}^{\pi}\sin^{n-2}(\psi)d\psi=2\int_{0}^{\frac{\pi}{2}}\sin^{n-2}(\psi)d\psi$.
The same derivation works when the conditional is on $\mathbf{U}=\mathbf{u}$:

$\displaystyle\mathbb{P}_{\mathbf{U}^{*}}\left\\{{\|\mathbf{U}^{*}{}^{\top}\mathbf{u}\|}\big{/}{(\|\mathbf{U}^{*}\|\|\mathbf{u}\|)}>\sqrt{1-\gamma_{2}^{2}}\bigg{|}\mathbf{u}\right\\}=\mathbb{P}_{\mathbf{U}^{*}}\left\\{|\cos(\psi)|>\sqrt{1-\gamma_{2}^{2}}\bigg{|}\mathbf{u}\right\\}$
(37)
$\displaystyle\qquad=\frac{2}{c_{1}}\int_{0}^{\arcsin\gamma_{2}}\sin^{n-2}(s)ds>\frac{\gamma_{2}^{n-2}\arcsin\gamma_{2}}{n-1},$

Because (8) and (37) do not involve $\mathbf{u}^{*}$ or $\mathbf{u}$, we have

$\displaystyle\mathbb{P}_{(\mathbf{U}^{*},\mathbf{U})}\big{\\{}\rho(\mathbf{U}^{*},\mathbf{U})>1-\gamma_{2}^{2}\big{\\}}$
$\displaystyle=\mathbb{P}_{\epsilon}\left\\{\rho(\mathbf{U}^{*},\mathbf{U})>1-\gamma_{2}^{2}\bigg{|}\mathbf{U}\right\\}$
$\displaystyle=\mathbb{P}_{\mathbf{U}}\left\\{\rho(\mathbf{U}^{*},\mathbf{U})>1-\gamma_{2}^{2}\bigg{|}\epsilon\right\\}>\frac{\gamma_{2}^{n-2}\arcsin\gamma_{2}}{n-1}.$

The above statement also suggests that $\rho(\mathbf{U}^{*},\mathbf{U})$ and
$\mathbf{U}^{*}$ are independent, therefore $\rho(\mathbf{U}^{*},\mathbf{U})$
and $\rho_{\tau^{\bot}}(\mathbf{U}^{,}\mathbf{X}_{0}\bm{\beta}_{0})$ are
independent.

∎

######  Lemma 9.

Let

$E(\gamma_{1},\gamma_{2})=\left\\{\max_{\tau\neq\tau_{0},|\tau|\leq|\tau_{0}|}\rho_{\tau^{\bot}}(\mathbf{U}^{*},\mathbf{X}_{0}\bm{\beta}_{0})<\gamma_{1}^{2},\rho(\mathbf{U}^{*},\mathbf{U})>1-\gamma_{2}^{2}\right\\},$

for any $0<\gamma_{1},\gamma_{2}<1$ and
$\rho(\mathbf{u},\tau)=\frac{\|\mathbf{H}_{\tau}\mathbf{u}\|^{2}}{\|\mathbf{u}\|^{2}}.$
Then

$\displaystyle\left\\{\rho(\mathbf{U}^{*},\mathbf{U})>1-\gamma_{2}^{2},\max_{\tau\neq\tau_{0},|\tau|\leq|\tau_{0}|}\rho_{\tau^{\bot}}(\mathbf{U},\mathbf{X}_{0}\bm{\beta}_{0})<\tilde{\gamma}_{1}^{2},\max_{\tau\neq\tau_{0},|\tau|\leq|\tau_{0}|}\rho(\mathbf{U},\tau)<1-\gamma_{2}\right\\}$
$\displaystyle\subset E(\gamma_{1},\gamma_{2}),$

where
$\tilde{\gamma}_{1}=(1-\sqrt{\gamma_{2}})\gamma_{1}-\sqrt{2-2\sqrt{1-\gamma^{2}_{2}}}.$

###### Proof of Lemma 9.

Denote by
$g_{\tau}(\mathbf{U})=\frac{\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}\|}{\|\mathbf{U}\|}$
and
$g_{\tau}(\mathbf{U}^{*})=\frac{\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}^{*}\|}{\|\mathbf{U}^{*}\|}$.
For each $\mathbf{U}^{*}$, given
$\left\\{\mathbf{U}\in\bigcap_{\\{\tau:|\tau|\leq|\tau_{0}|\\}}A(\tilde{\gamma}_{1}^{2},\tau)\right\\}$,
we have

$\displaystyle\frac{1}{\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}^{*}\|}((\mathbf{U}^{*})^{\top}(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{X}_{0}\bm{\beta}_{0}$
$\displaystyle=\frac{1}{\|\mathbf{U}\|}\mathbf{U}^{\top}(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{X}_{0}\bm{\beta}_{0}+\left((\mathbf{U}^{*})^{\top}/\|\mathbf{U}^{*}\|-\mathbf{U}^{\top}/\|\mathbf{U}\|\right)(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{X}_{0}\bm{\beta}_{0}$
$\displaystyle\quad+\left(\frac{1}{\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}^{*}\|}-\frac{1}{\|\mathbf{U}^{*}\|}\right)(\mathbf{U}^{*})^{\top}(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{X}_{0}\bm{\beta}_{0}$
$\displaystyle\leq\frac{\|\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}\|}{\|\mathbf{U}\|}\frac{1}{\|(\mathbf{I}-\mathbf{H}\tau)\mathbf{U}\|}\mathbf{U}^{\top}(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{X}_{0}\bm{\beta}_{0}+\left\|\frac{(\mathbf{U}^{*})^{\top}}{\|\mathbf{U}^{*}\|}-\frac{\mathbf{U}^{\top}}{\|\mathbf{U}\|}\right\|\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{X}_{0}\bm{\beta}_{0}\|$
$\displaystyle\quad+\frac{\|\mathbf{U}^{*}\|-\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}^{*}\|}{\|\mathbf{U}^{*}\|}\frac{1}{\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}^{*}\|}(\mathbf{U}^{*})^{\top}(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{X}_{0}\bm{\beta}_{0}$
$\displaystyle\leq
g_{\tau}(\mathbf{U})\tilde{\gamma}_{1}\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{X}_{0}\bm{\beta}_{0}\|+\sqrt{2-2\frac{(\mathbf{U}^{*})^{\top}\mathbf{U}}{\|(\mathbf{U}^{*})^{\top}\|\|\mathbf{U}\|}}\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{X}_{0}\bm{\beta}_{0}\|$
$\displaystyle\quad+(1-g_{\tau}(\mathbf{U}^{*}))\frac{1}{\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}^{*}\|}(\mathbf{U}^{*})^{\top}(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{X}_{0}\bm{\beta}_{0}$
$\displaystyle\leq
g_{\tau}(\mathbf{U})\tilde{\gamma}_{1}\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{X}_{0}\bm{\beta}_{0}\|+\sqrt{2-2\sqrt{1-\gamma^{2}_{2}}}\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{X}_{0}\bm{\beta}_{0}\|$
$\displaystyle\quad+(1-g_{\tau}(\mathbf{U}^{*}))\frac{1}{\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}^{*}\|}(\mathbf{U}^{*})^{\top}(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{X}_{0}\bm{\beta}_{0}.$

It then follows that

$\displaystyle\frac{1}{\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}^{*}\|}(\mathbf{U}^{*})^{\top}(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{X}_{0}\bm{\beta}_{0}$
$\displaystyle\leq\frac{g_{\tau}(\mathbf{U})}{g_{\tau}(\mathbf{U}^{*})}\tilde{\gamma}_{1}\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{X}_{0}\bm{\beta}_{0}\|+\frac{1}{g_{\tau}(\mathbf{U}^{*})}\sqrt{2-2\sqrt{1-\gamma^{2}_{2}}}\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{X}_{0}\bm{\beta}_{0}\|$

Further because
$\|(\mathbf{I}-\mathbf{H}\tau)\mathbf{U}^{*}\|\leq\|(\mathbf{I}-\mathbf{H}\tau\mathbb{P}_{\mathbf{U}})\mathbf{U}^{*}\|$,
if $(\mathbf{U}^{*},\mathbf{U})\in B(\gamma_{2})$ we have

$\displaystyle g_{\tau}(\mathbf{U}^{*})$
$\displaystyle\leq\frac{\|(I-\mathbb{P}_{\mathbf{U}})\mathbf{U}^{*}\|}{\|\mathbf{U}^{*}\|}+\frac{\|(\mathbb{P}_{\mathbf{U}}-\mathbf{H}_{\tau}\mathbb{P}_{\mathbf{U}})\mathbf{U}^{*}\|}{\|\mathbf{U}^{*}\|}$
$\displaystyle\leq\gamma_{2}+\frac{\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbb{P}_{\mathbf{U}}\mathbf{U}^{*}\|}{\|\mathbb{P}_{\mathbf{U}}\mathbf{U}^{*}\|}=\gamma_{2}+g_{\tau}(\mathbf{U}).$

Similarly, we can show that $g_{\tau}(\mathbf{U})\leq
g_{\tau}(\mathbf{U}^{*})+\gamma_{2}.$ Then

$\displaystyle\frac{g_{\tau}(\mathbf{U}^{*})}{g_{\tau}(\mathbf{U})}\geq
1-\frac{\gamma_{2}}{g_{\tau}(\mathbf{U})}.$

It then follows that a sufficient condition for
$\frac{1}{\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{U}^{*}\|}(\mathbf{U}^{*})^{\top}(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{X}_{0}\bm{\beta}_{0}\leq\gamma_{1}$
is

$\displaystyle\tilde{\gamma}_{1}\leq\left(1-\frac{\gamma_{2}}{g_{\tau}(\mathbf{U})}\right)\gamma_{1}-\sqrt{2-2\sqrt{1-\gamma^{2}_{2}}}\leq\frac{g_{\tau}(\mathbf{U}^{*})}{g_{\tau}(\mathbf{U})}\gamma_{1}-\sqrt{2-2\sqrt{1-\gamma^{2}_{2}}}$

Because it follows from $\rho(\mathbf{U},\tau)<1-\gamma_{2}$ that
$g_{\tau}(\mathbf{U})=\sqrt{1-\rho^{2}(\mathbf{U},\tau)}>\sqrt{\gamma_{2}},$
the above holds for
$\tilde{\gamma}_{1}=(1-\sqrt{\gamma_{2}})\gamma_{1}-\sqrt{2-2\sqrt{1-\gamma^{2}_{2}}}.$
Therefore

$\displaystyle\\{\rho(\mathbf{U}^{*},\mathbf{U})>1-\gamma_{2}^{2}\\}\bigcap\bigg{\\{}\max_{\tau\neq\tau_{0},|\tau|\leq|\tau_{0}|}\rho_{\tau^{\bot}}(\mathbf{U}^{*},\mathbf{X}_{0}\bm{\beta}_{0})<\tilde{\gamma}_{1}^{2},\max_{\tau\neq\tau_{0},|\tau|\leq|\tau_{0}|}\rho(\mathbf{U},\tau)<1-\gamma_{2}\bigg{\\}}$
$\displaystyle\subset E(\gamma_{1},\gamma_{2}).$

∎

######  Lemma 10.

Suppose $n-|\tau_{0}|>4.$ For any $0<\tilde{\gamma}_{1},\gamma_{2}<1,$

$\displaystyle\mathbb{P}\left(\mathbf{U}\not\in\left\\{\max_{\tau\neq\tau_{0},|\tau|\leq|\tau_{0}|}\rho_{\tau^{\bot}}(\mathbf{U},\mathbf{X}_{0}\bm{\beta}_{0})<\tilde{\gamma}_{1}^{2},\max_{\tau\neq\tau_{0},|\tau|\leq|\tau_{0}|}\rho(\mathbf{U},\tau)<1-\gamma_{2}\right\\}\right)$
$\displaystyle\leq\mathbb{P}\left(\max_{\tau\neq\tau_{0},|\tau|\leq|\tau_{0}|}\rho_{\tau^{\bot}}(\mathbf{U},\mathbf{X}_{0}\bm{\beta}_{0})\geq\tilde{\gamma}_{1}^{2}\right)+\mathbb{P}\left(\max_{\tau\neq\tau_{0},|\tau|\leq|\tau_{0}|}\rho(\mathbf{U},\tau)\geq
1-\gamma_{2}\right)$ $\displaystyle\leq
4(\arccos\tilde{\gamma}_{1})^{n-|\tau_{0}|-1}p^{|\tau_{0}|}+2\gamma_{2}^{\frac{n-|\tau_{0}|}{2}-1}(\sqrt{n}p)^{|\tau_{0}|}.$

###### Proof of Lemma 10.

Let
$g_{\tau}^{2}(\mathbf{u})=\frac{\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{u}\|^{2}}{\|\mathbf{u}\|^{2}}=\frac{\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{u}\|^{2}}{\|(\mathbf{I}-\mathbf{H}_{\tau})\mathbf{u}\|^{2}+\|\mathbf{H}_{\tau}\mathbf{u}\|^{2}},$