text:
\section{Introduction}

Consider the linear regression model,
\begin{equation}\label{eq:linearmodel}
\bfy=\mathbf{X}\bm{\beta}+ \bm{\varepsilon},
\end{equation}
where $\bfy$ is an $n$-vector, $\mathbf{X}$ an $n\times p$ design matrix, $\bm{\beta}=(\beta_j)_{1:p}$ the vector
of coefficients, and $\bm{\varepsilon}$ i.i.d. random errors with mean zero and variance $\sigma^2$.
Recently, $\ell_1$-penalized estimation methods \citep{Tibshirani96, Chen98}
have been widely used to find sparse estimates of the coefficient vector.
Given positive weights $w_j$, $j=1,\ldots,p$, and a tuning parameter $\lambda>0$,
an $\ell_1$-penalized estimator $\hat{\bmbeta}=(\hat{\beta}_j)_{1:p}$ is defined by minimizing
the following penalized loss function,
\begin{equation}\label{eq:lassoloss}
\ell(\bm{\beta})=\frac{1}{2}\|\bfy-\mathbf{X}\bm{\beta} \|_2^2 + n \lambda \sum_{j=1}^p w_j |\beta_j|.
\end{equation}
By different ways of choosing $w_j$, the estimator corresponds to the Lasso \citep{Tibshirani96}, the adaptive Lasso \citep{Zou06}, and the one-step linear local approximation (LLA) estimator \citep{Zou08} among others.
We call such an estimator a Lasso-type estimator.

In many applications of $\ell_1$-penalized regression, it is desired to quantify the uncertainty in the estimates.
However, except for very special cases, the sampling distribution of a Lasso-type estimator is complicated and
difficult to approximate. Closed-form approximations to the covariance matrices of
the estimators in \cite{Tibshirani96}, \cite{FanLi01}, and \cite{Zou06} are unsatisfactory,
as they all give zero variance for a zero component of the estimators and thus fail to quantify
the uncertainty in variable selection. Theoretical results on finite-sample distributions
and confidence sets of some
Lasso-type estimators have been developed \citep{Potscher09,Potscher10}
but only under orthogonal designs, which clearly limits general applications of these results.
The bootstrap can be used to approximate the sampling distribution of a Lasso-type estimator,
in which numerical optimization is needed to minimize \eqref{eq:lassoloss} for every bootstrap sample.
Although there are efficient algorithms, such as the Lars \citep{Efron04}, the homotopy algorithm \citep{Osborne00},
and coordinate descent \citep{Friedman07,WuLange08},
to solve this optimization problem, it is still time-consuming to apply these algorithms
hundreds or even thousands of times in bootstrap sampling.
As pointed out by \cite{KnightFu00} and \cite{Chatterjee10}, the bootstrap may not be consistent for estimating
the sampling distribution of the Lasso under certain circumstances.
To overcome this difficulty, a modified bootstrap \citep{Chatterjee11} and a perturbation resampling approach \citep{Minnier11} have been proposed, both justified under a fixed-$p$ asymptotic framework.
\cite{ZhangZhang11} have developed methods for constructing confidence intervals for individual coefficients and their linear combinations in high-dimensional ($p>n$) regression with sufficient conditions for the asymptotic normality of the proposed estimators. There are several recent articles on significance test and confidence region construction for sparse high-dimensional linear models \citep{Javanmard13b,Javanmard13a,Lockhart13,vandeGeer13}, all based on asymptotic distributions for various functions of the Lasso. On the other hand, knowledge on sampling distributions is also useful for distribution-based model selection with $\ell_1$ penalization, as demonstrated by stability selection \citep{Meinshausen10}
and the Bolasso \citep{Bach08}.

A possible alternative to the bootstrap or resampling
is to simulate from a sampling distribution
by Monte Carlo methods, such as Markov chain Monte Carlo (MCMC).
An obvious obstacle to using these methods for a Lasso-type estimator
is that its sampling distribution does not have a closed-form density. In this article, we study
the joint distribution of a Lasso-type estimator $\hat{\bmbeta}$ and the subgradient $\mathbf{S}$ of $\|\bm{\beta}\|_1$
evaluated at $\hat{\bmbeta}$.
Interestingly, this joint distribution has a density that can be calculated explicitly
assuming a normal error distribution, regardless of the relative size between $n$ and $p$.
Thus, one can develop Monte Carlo algorithms to draw samples from this joint distribution and estimate
various expectations of interest with respect to the sampling distribution of $\hat{\bmbeta}$, which is simply a marginal distribution. This approach offers great flexibility in studying the sampling distribution of
a Lasso-type estimator. For instance,
one may use importance sampling (IS) to accurately estimate a tail probability (small p-value) with respect
to the sampling distribution under a null hypothesis,
which can be orders of magnitude more efficient than any method directly targeting at the sampling distribution.
Another potential advantage of this approach is that, at each iteration, an MCMC algorithm only
evaluates a closed-form density, which is much faster than minimizing \eqref{eq:lassoloss}
numerically as used in the bootstrap. Furthermore, our method can be interpreted as an MCMC algorithm targeting at a multivariate normal distribution with locally reparameterized moves and hence is expected to be computationally tractable.

The remaining part of this article is organized as follows. After a high-level description of the basic idea, Section~\ref{sec:distribution} derives the density of the joint distribution of $\hat{\bmbeta}$ and $\mathbf{S}$ in the low-dimensional setting with $p\leq n$, and Section~\ref{sec:MC} develops MCMC algorithms for this setting. The density
in the high-dimensional setting with $p>n$ is derived in Section~\ref{sec:highdim}. In Section~\ref{sec:pvalue},
we construct applications of the high-dimensional result in p-value calculation
for Lasso-type inference by IS. Numerical examples are provided in Sections~\ref{sec:MC}
and \ref{sec:pvalue} to demonstrate the efficiency of the Monte Carlo algorithms. Section~\ref{sec:errorbounds} provides theoretical justifications for simulation from an estimated sampling distribution of the Lasso by establishing its  consistency as $p\gg n\to \infty$.
Section~\ref{sec:generalization} includes generalizations to random designs, a connection to model selection consistency,
and a Bayesian interpretation of the sampling distribution. The article concludes with a brief discussion and some remarks.
Technical proofs are relegated to Section~\ref{sec:proofs}.

Notations for vectors and matrices are defined here.
All vectors are regarded as column vectors.
Let $A=\{j_1,\ldots,j_k\}\subseteq \{1,\ldots,m\}$ and $B=\{i_1,\ldots,i_{\ell}\} \subseteq \{1,\ldots,n\}$ be two index sets.
For vectors $\mathbf{v}=(v_j)_{1:m}$ and $\mathbf{u}=(u_i)_{1:n}$,
we define $\mathbf{v}_A=(v_j)_{j\in A}=(v_{j_1},\ldots,v_{j_k})$, $\mathbf{v}_{-A}=(v_j)_{j\notin A}$,
and $(\mathbf{v}_A,\mathbf{u}_B)=(v_{j_1},\ldots,v_{j_k},u_{i_1},\ldots,u_{i_{\ell}})$.
For a matrix $\mathbf{M}=(M_{ij})_{m\times n}$,
write its columns as $\mathbf{M}_j$, $j=1,\ldots,n$. Then $\mathbf{M}_B=(\mathbf{M}_{j})_{j\in B}$ extracts the columns in $B$,
the submatrix $\mathbf{M}_{AB}=(M_{ij})_{i\in A, j \in B}$ extracts the rows in $A$ and the columns in $B$,
and $\mathbf{M}_{A\bullet}=(M_{ij})_{i \in A}$ extracts the rows in $A$. Furthermore, $\mathbf{M}^{\mathsf{T}}_B$ and
$\mathbf{M}_{AB}^{\mathsf{T}}$ are understood as $(\mathbf{M}_B)^{\mathsf{T}}$ and $(\mathbf{M}_{AB})^{\mathsf{T}}$, respectively.
We denote the row space, the null space, and the rank of $\mathbf{M}$ by $\text{row}(\mathbf{M})$, $\text{null}(\mathbf{M})$, and $\text{rank}(\mathbf{M})$, respectively.
Denote by $\mbox{diag}(\mathbf{v})$ the $m\times m$ diagonal matrix with $\mathbf{v}$ as the diagonal elements,
and by $\mbox{diag}(\mathbf{M},\mathbf{M}')$ the block diagonal matrix with $\mathbf{M}$ and $\mathbf{M}'$ as the diagonal blocks,
where the submatrices $\mathbf{M}$ and $\mathbf{M}'$ may be of different sizes and may not be square.
For a square matrix $\mathbf{M}$, $\mbox{diag}(\mathbf{M})$ extracts the diagonal elements.
Denote by $\mathbf{I}_n$ the $n\times n$ identity matrix.

\section{Estimator augmentation}\label{sec:distribution}

\subsection{The basic idea}\label{sec:basic}

Let $\mathbf{W}=\mbox{diag}(w_1,\ldots,w_p)$.
A minimizer $\hat{\bmbeta}$ of \eqref{eq:lassoloss} is given by the Karush-Kuhn-Tucker (KKT) condition
\begin{equation}\label{eq:lassograd}
\frac{1}{n}\mathbf{X}^{\mathsf{T}}\bfy = \frac{1}{n}\mathbf{X}^{\mathsf{T}} \mathbf{X} \hat{\bmbeta} + \lambda \mathbf{W} \mathbf{S},
\end{equation}
where $\mathbf{S}=(S_j)_{1:p}$ is the subgradient of the function $g(\bm{\beta})=\|\bm{\beta}\|_1$ evaluated at
the solution $\hat{\bmbeta}$. Therefore,
\begin{equation}\label{eq:defS}
\left\{
\begin{array}{ll}
S_j=\mbox{sgn}(\hat{\beta}_j) & \mbox{ if } \hat{\beta}_j \ne 0, \\
S_j \in [-1,1] & \mbox{ if } \hat{\beta}_j=0,
\end{array}
\right.
\end{equation}
for $j=1,\ldots,p$. Hereafter, we may simply call $\mathbf{S}$ the subgradient if the meaning is clear from context.
Lemma~\ref{lm:prelim} reviews a few basic facts about the uniqueness of $\hat{\bmbeta}$ and $\mathbf{S}$.
\begin{lemma}\label{lm:prelim}
For any $\bfy$, $\mathbf{X}$ and $\lambda>0$, every minimizer $\hat{\bmbeta}$ of \eqref{eq:lassoloss}
gives the same fitted value $\mathbf{X} \hat{\bmbeta}$ and the same subgradient $\mathbf{S}$. Moreover, if the columns of $\mathbf{X}$ are in general position, then $\hat{\bmbeta}$ is unique for any $\bfy$ and $\lambda>0$.
\end{lemma}
\begin{proof}
See Lemma 1 and Lemma 3 in \cite{Tibshirani13} for proof of the uniqueness of the fitted value $\mathbf{X}\hat{\bmbeta}$ and
the uniqueness of $\hat{\bmbeta}$. Since $\mathbf{S}$ is
a (vector-valued) function of $\mathbf{X}\hat{\bmbeta}$ from the KKT condition \eqref{eq:lassograd},
it is also unique for fixed $\bfy$, $\mathbf{X}$ and $\lambda$.
\end{proof}

We regard $\hat{\bmbeta}$ and $\mathbf{S}$ together as the solution to Equation \eqref{eq:lassograd}.
Lemma~\ref{lm:prelim} establishes that $(\hat{\bmbeta},\mathbf{S})$ is unique for any $\bfy$
assuming the columns of $\mathbf{X}$ are in general position (for a technical definition see \citealp{Tibshirani13}), regardless of the sizes of $n$ and $p$.
We call the vector $(\hat{\bmbeta},\mathbf{S})$ the augmented estimator in an $\ell_1$-penalized regression
problem. The augmented estimator will play a central role in our study of the sampling distribution of $\hat{\bmbeta}$.

Let $\mathbf{U} = \frac{1}{n}\mathbf{X}^{\mathsf{T}}\bm{\varepsilon}=\frac{1}{n}\mathbf{X}^{\mathsf{T}}\bfy-\mathbf{C}\bm{\beta}$, where
$\mathbf{C}=\frac{1}{n}\mathbf{X}^{\mathsf{T}} \mathbf{X}$ is the Gram matrix. By definition, $\mathbf{U}\in\text{row}(\mathbf{X})$. Rewrite the KKT condition as
\begin{equation}\label{eq:lassoKKTinU}
\mathbf{U} = \mathbf{C} \hat{\bmbeta} + \lambda \mathbf{W} \mathbf{S} - \mathbf{C} \bm{\beta} \defi \mathbf{H}(\hat{\bmbeta},\mathbf{S};\bm{\beta}),
\end{equation}
which shows that $\mathbf{U}$ is a function of $(\hat{\bmbeta},\mathbf{S})$. On the other hand,
$\bfy$ determines $(\hat{\bmbeta},\mathbf{S})$ only through $\mathbf{U}$, which implies that
$(\hat{\bmbeta},\mathbf{S})$ is unique for any $\mathbf{U}$ as long as it is unique for any $\bfy$.
Therefore, under the assumptions for the uniqueness of $\hat{\bmbeta}$,
$\mathbf{H}$ is a bijection between $(\hat{\bmbeta},\mathbf{S})$ and $\mathbf{U}$.
For a fixed $\mathbf{X}$, the only source of randomness in the linear model \eqref{eq:linearmodel}
is the noise vector $\bm{\varepsilon}$, which determines the distribution of $\mathbf{U}$.
With the bijection between $\mathbf{U}$ and $(\hat{\bmbeta},\mathbf{S})$, one may derive the joint distribution of $(\hat{\bmbeta},\mathbf{S})$,
which has a closed-form density under a normal error distribution. Then we develop Monte Carlo algorithms
to sample from this joint distribution and obtain the sampling distribution of $\hat{\bmbeta}$.
This is the key idea of this article, which works for both the low-dimensional setting ($p\leq n$) and the
high-dimensional setting ($p>n$). Although the basic strategy is the same, the technical details
are slightly more complicated for the high-dimensional setting.
For the sake of understanding, we first focus on the low-dimensional case in the remaining part of Section~\ref{sec:distribution} and Section~\ref{sec:MC}, and then generalize the results to the high-dimensional setting
in Section~\ref{sec:highdim}.

Before going through all the technical details, we take a glimpse of the utility of this work in a couple concrete examples. Given a design matrix $\mathbf{X}$ and a value of $\lambda$, our method gives a closed-form joint density $\pi$ for the Lasso-type estimator $\hat{\bmbeta}$ and the subgradient $\mathbf{S}$ under assumed values of the true parameters (Theorems~\ref{thm:joint} and \ref{thm:jointhigh}). Targeting at this density, we have developed MCMC algorithms, such as the Lasso sampler in Section~\ref{sec:MHLasso}, to draw samples from the joint distribution of $(\hat{\bmbeta},\mathbf{S})$. Such MCMC samples allow for approximation of marginal distributions for a Lasso-type estimator. Figure~\ref{fig:diagnostic} demonstrates the results of the Lasso sampler applied on a simulated dataset with $n=500$, $p=100$ and a normal error distribution. The scatter plot in Figure~\ref{fig:diagnostic}(a) confirms that $\hat{\beta}_j$ indeed may have a positive probability to be exactly zero in its sampling distribution. Accordingly, the distribution of the subgradient $S_j$ in Figure~\ref{fig:diagnostic}(b) has a continuous density on $(-1,1)$ and two point masses on $\pm 1$. The fast mixing and low autocorrelation shown in the figure are surprisingly satisfactory for a simple MCMC algorithm in such a high-dimensional and complicated space ($\mathbb{R}^p \times 2^{\{1,\ldots,p\}}$, see \eqref{eq:space}). Exploiting the explicit form of the bijection $\mathbf{H}$, we achieve the goal of sampling from the joint distribution of $(\hat{\bmbeta},\mathbf{S})$ via an MCMC algorithm essentially targeting at a multivariate normal distribution (Section~\ref{sec:viewofMCMC}). This approach does not need numerical optimization in any step, which makes it highly efficient compared to bootstrap or resampling-based methods. Another huge potential of our method is its ability to estimate tail probabilities, such as small p-values in a significance test (Section~\ref{sec:pvalue}). Estimating tail probabilities is challenging for any simulation method. With a suitable proposal distribution, having an explicit density makes it possible to accurately estimate tail probabilities by importance weights. For example, our method can estimate a tail probability on the order of $10^{-20}$, with a coefficient of variation around 2, by simulating only 5,000 samples from a proposal distribution. This is absolutely impossible when bootstrapping the Lasso or simulating from the sampling distribution directly.

\begin{figure}[t]
\centering
  \begin{minipage}[b]{0.4\linewidth}
   \centering
\includegraphics[width=\linewidth,trim=0in 0.2in 0in 0.5in,clip]{n500p100_s23_pair1v50.eps} \\
(a)
\end{minipage}
\begin{minipage}[b]{0.4\linewidth}
   \centering
   \includegraphics[width=\linewidth,trim=0in 0.2in 0in 0.5in,clip]{n500p100_s23_histS50.eps} \\
   (b)
\end{minipage}   \\
  \begin{minipage}[b]{0.4\linewidth}
   \centering
   \includegraphics[width=\linewidth,trim=0in 0.2in 0in 0.5in,clip]{n500p100_s23_beta1mix.eps} \\
   (c)
\end{minipage}
\begin{minipage}[b]{0.4\linewidth}
   \centering
   \includegraphics[width=\linewidth,trim=0in 0.2in 0in 0.5in,clip]{n500p100_s23_acf.eps} \\
   (d)
\end{minipage}   \\
   \caption{Demonstration of the Lasso sampler on a dataset with $p=100$. (a) Scatter plot of the samples of $\hat{\beta}_1$ (x-axis) and $\hat{\beta}_{50}$ (y-axis); (b) histogram of the subgradient $S_{50}$; (c) two sample paths of $\hat{\beta}_1$ with diverse initial values; (d) a typical autocorrelation function.\label{fig:diagnostic}}
\end{figure}

\subsection{The bijection}

In the low-dimensional setting, we assume that $\text{rank}(\mathbf{X})=p\leq n$, which guarantees that the columns of $\mathbf{X}$ are in general position.

Before writing down the bijection explicitly, we first examine the respective spaces for $\mathbf{U}$
and $(\hat{\bmbeta},\mathbf{S})$. Under the assumption that $\text{rank}(\mathbf{X})=p$, the row space of $\mathbf{X}$ is simply $\mathbb{R}^p$,
which is the space for $\mathbf{U}$. Let $\mathcal{A}=\text{supp}(\hat{\bmbeta})\defi\{j:\hat{\beta}_j \ne 0\}$ be the active set of $\hat{\bmbeta}$ and
$\mathcal{I}=\{1,\ldots,p\}\setminus \mathcal{A}$ be the inactive set, i.e., the set of the zero components of $\hat{\bmbeta}$.
After removing the degeneracies among its components as given in \eqref{eq:defS},
the vector $(\hat{\bmbeta},\mathbf{S})$ can be equivalently represented by the triple $(\hat{\bmbeta}_{\mathcal{A}},\mathbf{S}_{\mathcal{I}},\mathcal{A})$.
They are equivalent because
from $(\hat{\bmbeta}_{\mathcal{A}},\mathbf{S}_{\mathcal{I}},\mathcal{A})$ one can unambiguously recover $(\hat{\bmbeta},\mathbf{S})$,
by setting $\hat{\bmbeta}_{\mathcal{I}}=\mathbf{0}$ and $\mathbf{S}_{\mathcal{A}}=\mbox{sgn}(\hat{\bmbeta}_{\mathcal{A}})$ \eqref{eq:defS}, and vice versa.
It is more convenient and transparent to work with this equivalent representation.
One sees immediately that $(\hat{\bmbeta}_{\mathcal{A}},\mathbf{S}_{\mathcal{I}},\mathcal{A})$ lies in
\begin{equation}\label{eq:space}
\Omega= \{(\mathbf{b}_A,\mathbf{s}_I,A): A\subseteq\{1,\ldots,p\}, \mathbf{b}_A\in (\mathbb{R}\setminus\{0\})^{|A|}, \mathbf{s}_I \in [-1,1]^{p-|A|}\},
\end{equation}
where $I=\{1,\ldots,p\}\setminus A$. Hereafter, we always understand $(\mathbf{b}_A,\mathbf{s}_I,A)$ as the equivalent
representation of $(\mathbf{b},\mathbf{s})=((b_j)_{1:p},(s_j)_{1:p})$ with $\text{supp}(\mathbf{b})=A$ and $\mathbf{s}_A=\mbox{sgn}(\mathbf{b}_A)$.
Clearly, $\Omega \subset \mathbb{R}^p\times 2^{\{1,\ldots,p\}}$, where $2^{\{1,\ldots,p\}}$
is the collection of all subsets of $\{1,\ldots,p\}$, and thus $(\hat{\bmbeta}_{\mathcal{A}},\mathbf{S}_{\mathcal{I}},\mathcal{A})$
lives in the product space of $\mathbb{R}^p$ and a finite discrete space.

Partition $\hat{\bmbeta}$ as $(\hat{\bmbeta}_\mathcal{A},\hat{\bmbeta}_\mathcal{I})=(\hat{\bmbeta}_\mathcal{A},\mathbf{0})$ and
$\mathbf{S}$ as $(\mathbf{S}_\mathcal{A},\mathbf{S}_\mathcal{I})=(\mbox{sgn}(\hat{\bmbeta}_\mathcal{A}),\mathbf{S}_\mathcal{I})$.
Then the KKT condition \eqref{eq:lassoKKTinU} can be rewritten,
\begin{eqnarray}
\mathbf{U} & = &(\mathbf{C}_\mathcal{A} \mid \mathbf{C}_{\mathcal{I}})
\left(\begin{array}{c}
\hat{\bmbeta}_{\mathcal{A}} \\
\mathbf{0} \\
\end{array}\right)
+ \lambda (\mathbf{W}_{\mathcal{A}} \mid \mathbf{W}_{\mathcal{I}})
\left(\begin{array}{c}
\mathbf{S}_{\mathcal{A}} \\
\mathbf{S}_{\mathcal{I}} \\
\end{array}\right)
-\mathbf{C}\bm{\beta}, \label{eq:lassojoint_a} \\
& = & \mathbf{D}(\mathcal{A})\left(
\begin{array}{l}
\hat{\bmbeta}_{\mathcal{A}}  \\
\mathbf{S}_{\mathcal{I}}
\end{array}
\right) +  \lambda \mathbf{W}_{\mathcal{A}} \mbox{sgn}(\hat{\bmbeta}_\mathcal{A})-  \mathbf{C}\bm{\beta} \defi \mathbf{H}(\hat{\bmbeta}_\mathcal{A},\mathbf{S}_\mathcal{I}, \mathcal{A}; \bm{\beta}),
\label{eq:lassojoint}
\end{eqnarray}
where $\mathbf{D}(\mathcal{A})=(\mathbf{C}_{\mathcal{A}} \mid \lambda\mathbf{W}_{\mathcal{I}})$ is a $p\times p$ matrix. Permuting the rows of $\mathbf{D}(\mathcal{A})$,
one sees that
\begin{equation}\label{eq:defD}
|\det\mathbf{D}(\mathcal{A})|=
\det\left(
\begin{array}{cc}
\mathbf{C}_{\mathcal{A}\calA} & \mathbf{0} \\
\mathbf{C}_{\mathcal{I}\mathcal{A}} & \lambda\mathbf{W}_{\mathcal{I}\calI}
\end{array}
\right) = \lambda^{|\mathcal{I}|}\det(\mathbf{C}_{\mathcal{A}\calA})\prod_{j\in \mathcal{I}}w_j>0
\end{equation}
if $\mathbf{C}_{\mathcal{A}\calA}>0$.
Due to the equivalence between $(\hat{\bmbeta}_\mathcal{A},\mathbf{S}_\mathcal{I}, \mathcal{A})$ and $(\hat{\bmbeta},\mathbf{S})$,
the map $\mathbf{H}$ defined here is essentially
the same as the one defined in \eqref{eq:lassoKKTinU}.
\begin{lemma}\label{lm:bijectioninlowdim}
If $\text{rank}(\mathbf{X})=p$, then for any $\bm{\beta}$ and $\lambda>0$, the mapping $\mathbf{H}:\Omega\to \mathbb{R}^p$ defined in \eqref{eq:lassojoint} is a bijection that maps $\Omega$ onto $\mathbb{R}^p$.
\end{lemma}
\begin{proof}
For any $\mathbf{U} \in \mathbb{R}^p$, there is a unique solution $(\hat{\bmbeta},\mathbf{S})$ to Equation \eqref{eq:lassoKKTinU} if $\text{rank}(\mathbf{X})=p$,
and thus, a unique $(\hat{\bmbeta}_\mathcal{A},\mathbf{S}_\mathcal{I},\mathcal{A})\in\Omega$
such that $\mathbf{H}(\hat{\bmbeta}_\mathcal{A},\mathbf{S}_\mathcal{I}, \mathcal{A}; \bm{\beta})=\mathbf{U}$.
For any $(\hat{\bmbeta}_\mathcal{A},\mathbf{S}_\mathcal{I},\mathcal{A}) \in \Omega$, $\mathbf{H}$ maps it into $\mathbb{R}^p$.
\end{proof}

It is helpful for understanding the map $\mathbf{H}$ to consider its inverse $\mathbf{H}^{-1}$ and its restriction to $\mathcal{A}=A$,
where $A$ is a fixed subset of $\{1,\ldots,p\}$.
For any $\mathbf{U} \in \mathbb{R}^p$,
if $\mathbf{H}^{-1}(\mathbf{U}; \bm{\beta}) = (\hat{\bmbeta}_\mathcal{A},\mathbf{S}_\mathcal{I}, \mathcal{A})$, then the unique solution to Equation~\eqref{eq:lassojoint_a} is $(\hat{\bmbeta}_\mathcal{A},\mathbf{S}_\mathcal{I}, \mathcal{A})$.
Given a fixed $A$, $(\hat{\bmbeta}_A,\mathbf{S}_I)$ lives in the subspace
\begin{equation}\label{eq:subspace}
\Omega_A=\{(\mathbf{b}_A,\mathbf{s}_I) \in \mathbb{R}^p: \mathbf{b}_A\in (\mathbb{R}\setminus\{0\})^{|A|}, \mathbf{s}_I \in [-1,1]^{p-|A|}\}.
\end{equation}
Let $\mathbf{H}_A(\mathbf{b}_A,\mathbf{s}_I;\bm{\beta}) = \mathbf{H}(\mathbf{b}_A,\mathbf{s}_I,A;\bm{\beta})$ for $(\mathbf{b}_A,\mathbf{s}_I)\in \Omega_A$
and $U_A=\mathbf{H}_A(\Omega_A;\bm{\beta})$
be the image of $\Omega_A$ under the map $\mathbf{H}_A$. Now imagine we plug different $\mathbf{U} \in \mathbb{R}^p$
into Equation~\eqref{eq:lassojoint_a} and solve for $(\hat{\bmbeta}_\mathcal{A},\mathbf{S}_\mathcal{I}, \mathcal{A})$.
Then the set $\Omega_A\times \{A\}$ is the collection of all
possible solutions such that $\text{supp}(\hat{\bmbeta})=A$, the set $U_A$ is the collection of all $\mathbf{U}$ that give
these solutions, and $\mathbf{H}_A$ is a bijection between the two sets. It is easy to see that
$\Omega=\bigcup_A \Omega_A \times \{A\}$,
i.e., $\{\Omega_A\times\{A\}\}$, for $A$ extending over all subsets of $\{1,\ldots,p\}$,
form a partition of the space $\Omega$. The bijective
nature of $\mathbf{H}$ implies that $\{U_A\}$ also form a partition of $\mathbb{R}^p$, the space of $\mathbf{U}$.
Figure~\ref{fig:bijection} illustrates the bijection $\mathbf{H}$ for $p=2$ and the space partitioning by $A$.
In this case, $\mathbf{H}_A$ map the four subspaces $\Omega_A$ for $A=\varnothing,\{1\},\{2\},\{1,2\}$,
each in a different $\mathbb{R}^2$, onto the space of $\mathbf{U}$ which is another $\mathbb{R}^2$.

\begin{figure}[ht]
\centering
   \includegraphics[width=0.55\linewidth,angle=-90,trim=0.5in 0.5in 0.5in 0.5in,clip]{Bijection.eps} \\
   \caption{The bijection $\mathbf{H}$, its restrictions $\mathbf{H}_A$, the four subspaces $\Omega_A$ (shaded areas) and
   the corresponding partition in the space of $\mathbf{U}$ for $p=2$.
    \label{fig:bijection}}
\end{figure}

\begin{remark}\label{rm:bijectionlowdim}
The simple fact that $\mathbf{H}$ maps every point in $\Omega$ into $\text{row}(\mathbf{X})=\mathbb{R}^p$
is crucial to the derivation of the sampling distribution of
$(\hat{\bmbeta}_{\mathcal{A}},\mathbf{S}_{\mathcal{I}},\mathcal{A})$ in the low-dimensional setting. This means that every
$(\mathbf{b}_A,\mathbf{s}_I,A) \in \Omega$ is the solution to Equation \eqref{eq:lassojoint_a}
for $\mathbf{U}=\mathbf{u}=\mathbf{H}(\mathbf{b}_A,\mathbf{s}_I,A;\bm{\beta})$, and therefore one can simply find the probability density
of $(\hat{\bmbeta}_{\mathcal{A}},\mathbf{S}_{\mathcal{I}},\mathcal{A})$ at $(\mathbf{b}_A,\mathbf{s}_I,A)$ by the density of $\mathbf{U}$ at  $\mathbf{u}$.
This is not the case when $p>n$ (Section~\ref{sec:highdim}).
\end{remark}

\subsection{The sampling distribution}

Now we can use the bijection $\mathbf{H}$ to find the distribution of $(\hat{\bmbeta}_\mathcal{A},\mathbf{S}_\mathcal{I}, \mathcal{A})$ from the distribution of $\mathbf{U}$. Let $\xi_k$ denote $k$-dimensional Lebesgue measure.

\begin{theorem}\label{thm:joint}
Assume that $\text{rank}(\mathbf{X})=p$ and let $f_{\mathbf{U}}$ be the probability density of $\mathbf{U}$ with respect to $\xi_p$.
For $(\mathbf{b}_A,\mathbf{s}_I,A)\in\Omega$, the joint distribution of $(\hat{\bmbeta}_\mathcal{A},\mathbf{S}_\mathcal{I}, \mathcal{A})$ is given by
\begin{eqnarray}\label{eq:jointdensity}
 P(\hat{\bmbeta}_A\in d\mathbf{b}_A,\mathbf{S}_I \in d\mathbf{s}_I,\mathcal{A}=A)
&=&f_{\mathbf{U}}(\mathbf{H}(\mathbf{b}_A,\mathbf{s}_I,A;\bm{\beta}))|\det\mathbf{D}(A)| \xi_p(d\mathbf{b}_A d\mathbf{s}_I) \nonumber \\
& \defi & \pi(\mathbf{b}_A,\mathbf{s}_I,A) \xi_p(d\mathbf{b}_A d\mathbf{s}_I),
\end{eqnarray}
and the distribution of $(\hat{\bmbeta}_\mathcal{A}, \mathcal{A})$ is a marginal distribution given by
\begin{align}\label{eq:distrnbeta}
 P(  \hat{\bmbeta}_A \in d\mathbf{b}_A, \mathcal{A}=A)
 =\left[\int_{[-1,1]^{p-|A|}}\pi(\mathbf{b}_A,\mathbf{s}_I,A)\xi_{p-|A|}(d\mathbf{s}_I)\right]\xi_{|A|}(d\mathbf{b}_A).
\end{align}
\end{theorem}
\begin{proof}
Let $\mathbf{u}=\mathbf{H}(\mathbf{b}_A,\mathbf{s}_I,A;\bm{\beta})= \mathbf{H}_A(\mathbf{b}_A,\mathbf{s}_I;\bm{\beta})$.
From \eqref{eq:lassojoint} and \eqref{eq:subspace}, one sees that
for any fixed $A$, $b_j\ne 0$ for all $j\in A$ and $\mathbf{H}_A$ is differentiable. Differentiating $\mathbf{u}$ with respect to $(\mathbf{b}_A,\mathbf{s}_I)$,
\begin{equation*}
d\mathbf{u}=\frac{\partial \mathbf{H}_A}{\partial(\mathbf{b}_A,\mathbf{s}_I)^{\mathsf{T}}}\left(\begin{array}{l}
d\mathbf{b}_A  \\ d\mathbf{s}_I \end{array}\right)
=\mathbf{D}(A) \left(\begin{array}{l}
d\mathbf{b}_A  \\ d\mathbf{s}_I \end{array}\right)
\end{equation*}
and thus $\xi_p(d\mathbf{u})=|\det\mathbf{D}(A)|\xi_p(d\mathbf{b}_A d\mathbf{s}_I)$.  Since $\mathbf{H}$ and
$\mathbf{H}_A:\Omega_A \to U_A$ are bijections, a change of variable gives
\begin{align*}
& P(\hat{\bmbeta}_A\in d\mathbf{b}_A,\, \mathbf{S}_I \in d\mathbf{s}_I,\mathcal{A}=A)  =  P(\mathbf{U} \in d\mathbf{u}) \\
   &\quad\quad =  f_{\mathbf{U}}(\mathbf{H}(\mathbf{b}_A,\mathbf{s}_I,A;\bm{\beta})) |\det\mathbf{D}(A)|\xi_p(d\mathbf{b}_A d\mathbf{s}_I).
\end{align*}
Integrating \eqref{eq:jointdensity} over $\mathbf{s}_I \in [-1,1]^{p-|A|}$
gives \eqref{eq:distrnbeta}.
\end{proof}

\begin{remark}\label{rm:compare}
Equation~\eqref{eq:jointdensity} gives the joint distribution of $(\hat{\bmbeta}_\mathcal{A},\mathbf{S}_\mathcal{I}, \mathcal{A})$ and
effectively the joint distribution of $(\hat{\bmbeta},\mathbf{S})$. The density $\pi(\mathbf{b}_A,\mathbf{s}_I,A)$ is defined
with respect to the product of $\xi_p$ and counting measure on $2^{\{1,\ldots,p\}}$. Analogously, the sampling distribution
of $\hat{\bmbeta}$ is given by the distribution of $(\hat{\bmbeta}_\mathcal{A}, \mathcal{A})$ in \eqref{eq:distrnbeta}. To be rigorous, \eqref{eq:jointdensity} is derived by assuming that $(\mathbf{b}_A,\mathbf{s}_I)$ is an interior point of $\Omega_A$. Note that $(\mathbf{b}_A,\mathbf{s}_I)\in\Omega_A$ is not in the interior if and only if $|s_j|=1$ for some $j\in I$, and thus the Lebesgue measure of the union of these points is zero. Therefore, it will cause no problem at all to use $\pi$ as the density for all points in $\Omega$. The joint distribution of $(\hat{\bmbeta}_\mathcal{A},\mathbf{S}_\mathcal{I}, \mathcal{A})$ has at least two nice properties which make it much more tractable than the distribution of $\hat{\bmbeta}$. First, the density $\pi$ does not involve multidimensional integral and has a closed-form expression that can be calculated explicitly if $f_{\mathbf{U}}$ is given.
Second, the continuous components $(\hat{\bmbeta}_\mathcal{A},\mathbf{S}_\mathcal{I})$
always have the same dimension $(=p)$ for any value of $\mathcal{A}$,
while $\hat{\bmbeta}_\mathcal{A}$ lives in $\mathbb{R}^{|\mathcal{A}|}$ whose dimension changes
with $\mathcal{A}$. These two properties are critical to the development of MCMC to sample from $\pi$.
See Section~\ref{sec:reversibility} for more discussion. We explicitly include the dominating Lebesgue measure to clarify
the dimension of a density.
\end{remark}

\begin{remark}\label{rm:marginal}
The distribution of $\hat{\bmbeta}=(\hat{\beta}_1,\ldots,\hat{\beta}_p)$ in \eqref{eq:distrnbeta} is essentially
defined for each $A$. In many problems, one may be interested in the marginal distribution of $\hat{\beta}_j$
such as for calculating p-values and constructing confidence intervals. To obtain such a marginal distribution,
we need to sum over all possible active sets, which cannot be done analytically. Our strategy is to
draw samples from the joint distribution of $(\hat{\bmbeta}_\mathcal{A}, \mathbf{S}_\mathcal{I},\mathcal{A})$ by a Monte Carlo method. Then from
the Monte Carlo samples one can easily approximate any marginal distribution of interest, such as that of $\hat{\beta}_j$.
This is exactly our motivation
for estimator augmentation, which is in spirit similar to the use of auxiliary variables in the MCMC literature.
\end{remark}

To further help our understanding of the density $\pi$,
consider a few conditional and marginal distributions derived from the joint
distribution \eqref{eq:jointdensity}. First, the sampling distribution of the active set $\mathcal{A}$ is given by
\begin{equation}\label{eq:distrnforA}
P(\mathcal{A}=A)=\int_{\Omega_A} \pi(\mathbf{b}_A,\mathbf{s}_I,A) \xi_p(d\mathbf{b}_A d\mathbf{s}_I)\defi Z_A,
\end{equation}
where $\Omega_A$ is the subspace for $(\hat{\bmbeta}_A,\mathbf{S}_I)$ defined in \eqref{eq:subspace}. In other words,
$Z_A$ is the probability of $\Omega_A\times\{A\}$ with respect to the joint distribution $\pi$.
Second, the conditional density of $(\hat{\bmbeta}_A,\mathbf{S}_I)$ given $\mathcal{A}=A$ (with respect to $\xi_p$) is
\begin{equation}\label{eq:condforbs}
\pi(\mathbf{b}_A,\mathbf{s}_I \mid A)=\frac{1}{Z_A} \pi(\mathbf{b}_A,\mathbf{s}_I,A)\propto f_{\mathbf{U}}(\mathbf{H}(\mathbf{b}_A,\mathbf{s}_I,A;\bm{\beta}))
\end{equation}
for $(\mathbf{b}_A,\mathbf{s}_I)\in \Omega_A \subset \mathbb{R}^p$.
Using $p=2$ as an illustration, the joint density
$\pi$ is defined over all four shaded areas in Figure~\ref{fig:bijection},
while a conditional density $\pi(\cdot \mid A)$ is defined on each one of them.
To give a concrete probability calculation, for $a_2>a_1>0$,
\begin{eqnarray*}
P(\hat{\beta}_1 \in [a_1,a_2], \hat{\beta}_2=0) & = & P(\hat{\beta}_1 \in[a_1,a_2], \mathcal{A}=\{1\}) \\
& = & \int_{a_1}^{a_2} \int_{-1}^{1} \pi(b_1,s_2,\{1\}) ds_2 db_1,
\end{eqnarray*}
which is an integral over the rectangle $[a_1,a_2] \times [-1,1]$ in $\Omega_{\{1\}}$ (Figure~\ref{fig:bijection}).
Clearly, this probability can be approximated by Monte Carlo integration if we have enough samples from $\pi$.

\begin{remark}
We emphasize that the weights $w_j$ and the tuning parameter $\lambda$ are assumed to be fixed in Theorem~\ref{thm:joint}. For the adaptive Lasso, one may choose $w_j=|\tilde{\beta}_j|^{-1}$ based on an initial estimate $\tilde{\beta}_j$. The distribution \eqref{eq:jointdensity} is valid for the adaptive Lasso only if we ignore the randomness in $\tilde{\beta}_j$ and regard $w_j$ as constants during the repeated sampling procedure.
\end{remark}

\subsection{Normal errors}

Denote by $\mathcal{N}_k(\bm{\mu},\bm{\Sigma})$ the $k$-variate normal distribution with mean $\bm{\mu}$ and covariance
matrix $\bm{\Sigma}$, and by $\phi_k(\bfz;\bm{\mu},\bm{\Sigma})$ its probability density function.
If the error $\bm{\varepsilon} \sim \mathcal{N}_n(\mathbf{0},\sigma^2 \mathbf{I}_n)$ and $\text{rank}(\mathbf{X})=p$,
then $\mathbf{U} \sim \mathcal{N}_p(\mathbf{0},\frac{\sigma^2}{n}\mathbf{C})$.
In this case, the joint density $\pi$ \eqref{eq:jointdensity} has a closed-form expression.
Recall that $\mathbf{s}_A=\mbox{sgn}(\mathbf{b}_A)$ and define
\begin{eqnarray}
 \bm{\mu} (A,\mathbf{s}_A;\bm{\beta}) & = & [\mathbf{D}(A)]^{-1} \left( \mathbf{C}\bm{\beta} - \lambda  \mathbf{W}_A \mathbf{s}_A\right),  \label{eq:munorm} \\
\bm{\Sigma}(A;\sigma^2) & = & \frac{\sigma^2}{n}[\mathbf{D}(A)]^{-1}\mathbf{C}[\mathbf{D}(A)]^{-\mathsf{T}}.  \label{eq:varnorm}
\end{eqnarray}
\begin{corollary}\label{cor:normaljoint}
If $\text{rank}(\mathbf{X})=p$ and $\bm{\varepsilon} \sim \mathcal{N}_n(\mathbf{0},\sigma^2 \mathbf{I}_n)$, then
the joint density of $(\hat{\bmbeta}_{\mathcal{A}},\mathbf{S}_{\mathcal{I}},\mathcal{A})$ is
\begin{equation}\label{eq:jointdensitynormal}
\pi(\mathbf{b}_A,\mathbf{s}_I,A) =\phi_p(\bfz ; \bm{\mu}(A,\mathbf{s}_A;\bm{\beta}),\bm{\Sigma}(A;\sigma^2))\bm{1}((\bfz,A)\in \Omega),
\end{equation}
where $\bfz=(\mathbf{b}_A,\mathbf{s}_I)\in\mathbb{R}^p$ and $\bm{1}(\cdot)$ is an indicator function.
\end{corollary}
\begin{proof}
First note that
\begin{equation}\label{eq:affineKKT}
\mathbf{H}(\mathbf{b}_A,\mathbf{s}_I,A;\bm{\beta})=\mathbf{D}(A)[\bfz-\bm{\mu}(A,\mathbf{s}_A;\bm{\beta})].
\end{equation}
Under the assumptions, $\mathbf{U} \sim \mathcal{N}_p(\mathbf{0},\frac{\sigma^2}{n}\mathbf{C})$. By Theorem~\ref{thm:joint},
\begin{eqnarray*}
\pi(\mathbf{b}_A,\mathbf{s}_I,A)&=&\phi_p\left(\mathbf{D}(A)[\bfz-\bm{\mu}(A,\mathbf{s}_A;\bm{\beta})]; \mathbf{0},n^{-1}{\sigma^2}\mathbf{C}\right)|\det\mathbf{D}(A)| \\
& = & \phi_p\left(\bfz ; \bm{\mu}(A,\mathbf{s}_A;\bm{\beta}),n^{-1}{\sigma^2}[\mathbf{D}(A)]^{-1}\mathbf{C}[\mathbf{D}(A)]^{-\mathsf{T}}\right) \\
& = & \phi_p(\bfz ; \bm{\mu}(A,\mathbf{s}_A;\bm{\beta}),\bm{\Sigma}(A;\sigma^2))
\end{eqnarray*}
for $(\mathbf{b}_A,\mathbf{s}_I,A)=(\bfz,A)\in \Omega$.
\end{proof}

Without the normal error assumption,
Corollary~\ref{cor:normaljoint} is still a good approximation when $n$ is large and $p$ is fixed,
since $\sqrt{n} \mathbf{U} \toinL \mathcal{N}_p(\mathbf{0},\sigma^2\mathbf{C})$ assuming $\frac{1}{n}\mathbf{X}^{\mathsf{T}}\mathbf{X}\to \mathbf{C}>0$ as $n\to\infty$.

Note that both the continuous components $\bfz$ and the active set $A$ are
arguments of the density \eqref{eq:jointdensitynormal}.
For different $A$ and $\mathbf{s}_A$, the normal density $\phi_p$ has different parameters.
Given a particular $A^*$ and $\mathbf{s}^*\in \{\pm 1\}^{|A^*|}$, let $I^*=\{1,\ldots,p\}\setminus A^*$ and
\begin{equation}\label{eq:subsetAs}
\Omega_{A^*,\mathbf{s}^*}=\{(\mathbf{b}_{A^*},\mathbf{s}_{I^*})\in\Omega_{A^*}: \mbox{sgn}(\mathbf{b}_{A^*})=\mathbf{s}^*\}.
\end{equation}
Then $\Omega_{A^*,\mathbf{s}^*}\times\{A^*\}$ is the subset of $\Omega$ corresponding to
the event $\{\mathcal{A}=A^*, \mbox{sgn}(\hat{\bmbeta}_{A^*})=\mathbf{s}^*\}$.
For $\bfz \in\Omega_{A^*,\mathbf{s}^*}$, the density $\pi(\bfz,A^*)$ is identical to
$\phi_p(\bfz ; \bm{\mu}(A^*,\mathbf{s}^*;\bm{\beta}),\bm{\Sigma}(A^*;\sigma^2))$, i.e.,
\begin{equation*}
\pi(\bfz,A^*)\bm{1}(\bfz \in\Omega_{A^*,\mathbf{s}^*})
=\phi_p(\bfz ; \bm{\mu}(A^*,\mathbf{s}^*;\bm{\beta}),\bm{\Sigma}(A^*;\sigma^2))\bm{1}(\bfz \in\Omega_{A^*,\mathbf{s}^*}).
\end{equation*}
Intuitively, this is because $\mathbf{H}$ restricted to $A=A^*$ and $\mathbf{s}_{A^*}=\mathbf{s}^*$ is simply an affine map
[see \eqref{eq:affineKKT}].
Consequently, the probability of $\Omega_{A^*,\mathbf{s}^*}\times\{A^*\}$ with respect to $\pi$ is
\begin{equation}\label{eq:signpm}
P(\mathcal{A}=A^*, \mbox{sgn}(\hat{\bmbeta}_{A^*})=\mathbf{s}^*)=\int_{\Omega_{A^*,\mathbf{s}^*}} \phi_p(\bfz ; \bm{\mu}(A^*,\mathbf{s}^*;\bm{\beta}),\bm{\Sigma}(A^*;\sigma^2)) \xi_p(d\bfz),
\end{equation}
and $[\hat{\bmbeta}_{A^*},\mathbf{S}_{I^*} \mid \mathcal{A}=A^*,\mbox{sgn}(\hat{\bmbeta}_{A^*})=\mathbf{s}^*]$
is the truncated $\mathcal{N}_p(\bm{\mu}(A^*,\mathbf{s}^*;\bm{\beta}),\bm{\Sigma}(A^*;\sigma^2))$ on $\Omega_{A^*,\mathbf{s}^*}$.
For $p=2$, if $A^*=\{1\}$, and $\mathbf{s}^*=-1$, the region $\Omega_{\{1\},-1}=(-\infty,0)\times[-1,1]$ is
the left half of the $\Omega_{\{1\}}$ in Figure~\ref{fig:bijection} and the density $\pi$ restricted to this region is
the same as the part of a bivariate normal density on the same region.

If $\mathbf{C}=\mathbf{I}_p$ and $\mathbf{W}=\mathbf{I}_p$, the Lasso is equivalent to soft-thresholding the ordinary least-squares estimator $\hat{\bmbeta}^{\text{OLS}}=(\hat{\beta}^{\text{OLS}}_j)_{1:p}$. In this case, \eqref{eq:munorm} and \eqref{eq:varnorm} have simpler forms:
\begin{eqnarray*}
\bm{\mu} (A,\mathbf{s}_A;\bm{\beta})& =&
\left(
\begin{array}{c}
\bm{\beta}_{A}-\lambda \mathbf{s}_{A}  \\
\lambda^{-1} \bm{\beta}_I
\end{array}
\right), \\
\bm{\Sigma}(A;\sigma^2)& = &
\frac{\sigma^2}{n}\left( \begin{array}{cc}
\mathbf{I}_{|A|} &  \mathbf{0} \\
\mathbf{0} & \lambda^{-2} \mathbf{I}_{|I|}
 \end{array}\right).
\end{eqnarray*}
By \eqref{eq:signpm} we find, for example,
\begin{align*}
&P(\mathcal{A}=A,\mbox{sgn}(\hat{\bmbeta}_A)=(1,\ldots,1)) \\
& = \prod_{j\in A} \int_{0}^{\infty} \phi\left(b_j;\beta_j-\lambda,\frac{\sigma^2}{n}\right) db_j \cdot
\prod_{j\in I} \int_{-1}^1 \phi\left(s_j;\frac{\beta_j}{\lambda},\frac{\sigma^2}{\lambda^2 n}\right) ds_j \\
&= \prod_{j\in A} P(\hat{\beta}^{\text{OLS}}_j> \lambda) \cdot  \prod_{j\in I} P(|\hat{\beta}^{\text{OLS}}_j| \leq \lambda),
\end{align*}
where the last equality is due to that $\hat{\bmbeta}^{\text{OLS}} \sim \mathcal{N}_p(\bm{\beta},n^{-1}\sigma^2 \mathbf{I}_p)$. One sees that our result is consistent with that  obtained directly from soft-thresholding each component of $\hat{\bmbeta}^{\text{OLS}}$ by $\lambda$.

\subsection{Estimation}\label{sec:asymest}

To apply Theorem~\ref{thm:joint} in practice, one needs to estimate $f_{\mathbf{U}}$ and $\bm{\beta}$ if they are not given.
Suppose that $f_{\mathbf{U}}$ is estimated by $\hat{f}_{\mathbf{U}}$ and $\bm{\beta}$ is estimated by
$\check{\bmbeta}$. Then, the corresponding estimate of the density $\pi$ is
\begin{equation}\label{eq:estjointdensity}
\hat{\pi}(\mathbf{b}_A,\mathbf{s}_I,A) = \hat{f}_{\mathbf{U}}(\mathbf{H}(\mathbf{b}_A,\mathbf{s}_I,A;\check{\bmbeta}))|\det\mathbf{D}(A)|.
\end{equation}

Since $\mathbb{E}(\mathbf{U})=\mathbf{0}$ and $\text{Var}(\sqrt{n}\mathbf{U})=\sigma^2\mathbf{C}$, estimating $f_{\mathbf{U}}$
reduces to estimating $\sigma^2$ when $\bm{\varepsilon}$ is normally distributed or when the sample size $n$ is large.
A consistent estimator of $\sigma^2$ can be constructed given a consistent estimator of $\bm{\beta}$.
For example, when $p<n$ one may use
\begin{equation}\label{eq:estsigma2}
\hat{\sigma}^2 = \frac{\| \bfy-\mathbf{X} \check{\bmbeta}\|_2^2}{n-p},
\end{equation}
provided that $\check{\bmbeta}$ is consistent for $\bm{\beta}$.
If $\bm{\varepsilon}$ does not follow a normal distribution,
one can apply other parametric or nonparametric methods to estimate $f_{\mathbf{U}}$.
Here, we propose a bootstrap-based approach under the assumption that $\mathbf{U}$ is elliptically symmetric.
That is, $\tilde{\mathbf{U}}=\mathbf{C}^{-1/2}\mathbf{U}$ is spherically symmetric: For $\mathbf{v}_1,\mathbf{v}_2 \in \mathbb{R}^p$,
if $\|\mathbf{v}_1\|_2=\|\mathbf{v}_2\|_2$ then $f_{\tilde{\mathbf{U}}}(\mathbf{v}_1)=f_{\tilde{\mathbf{U}}}(\mathbf{v}_2)$,
where $f_{\tilde{\mathbf{U}}}$ is the density of $\tilde{\mathbf{U}}$.
Generate bootstrap samples, $\bm{\varepsilon}^{(i)}=(\varepsilon^{(i)}_1,\ldots,\varepsilon^{(i)}_n)$ for
$i=1,\ldots,K$, by resampling with replacement from
$\hat{\bm{\varepsilon}}=(\hat{\varepsilon}_1,\ldots,\hat{\varepsilon}_n)=\bfy-\mathbf{X} \check{\bmbeta}$,
and calculate $\tilde{\mathbf{U}}^{(i)}=\frac{1}{n}\mathbf{C}^{-1/2}\mathbf{X}^{\mathsf{T}}\bm{\varepsilon}^{(i)}$ for each $i$.
Given $0=h_0 < h_1 < \cdots < h_M<\infty$, let $K_m=|\{i:h_{m-1} \leq \|\tilde{\mathbf{U}}^{(i)}\|_2 < h_m\}|$ for $m=1,\ldots,M$.
The density of $\tilde{\mathbf{U}}$ is then estimated by
\begin{equation}
\hat{f}_{\tilde{\mathbf{U}}}(\mathbf{v}) \propto \sum_{m=1}^M \frac{K_m}{h_m^p - h_{m-1}^p}  \bm{1} (h_{m-1} \leq \|\mathbf{v}\|_2 < h_m)
\end{equation}
for $\|\mathbf{v}\|_2 \in [0,h_M)$. The density for $\|\mathbf{v}\|_2 \geq h_M$ can be estimated by linear extrapolation
of $\log\hat{f}_{\tilde{\mathbf{U}}}$. Finally, set $\hat{f}_{\mathbf{U}}(\mathbf{u}) = \hat{f}_{\tilde{\mathbf{U}}}(\mathbf{C}^{-1/2}\mathbf{u})(\det\mathbf{C})^{-1/2}$.

In general, estimating $f_{\mathbf{U}}$ is difficult when $p$ is large. One may have to assume some parametric
density for $\mathbf{U}$, which reduces the problem to the estimation of a few unknown parameters.
Besides normality, one may assume that $\mathbf{U}$ follows a multivariate $t$ distribution, which
is motivated from a Bayesian perspective to be discussed in Section~\ref{sec:Bayesian}.

Sampling from $\pi$ or $\hat{\pi}$ can be very useful for statistical inference based on a Lasso-type estimator. We may directly draw $({\bmbeta}^*,\mathbf{S}^*)$ from $\hat{\pi}$ given $(\check{\bmbeta},\hat{\sigma})$ and use the conditional distribution $[({\bmbeta}^*-\check{\bmbeta})\mid \check{\bmbeta},\hat{\sigma}]$ to construct confidence regions around $\hat{\bmbeta}$. Under some assumptions, the conditional distribution $[({\bmbeta}^*-\check{\bmbeta})\mid \check{\bmbeta},\hat{\sigma}]$ provides a valid approximation to the true sampling distribution of $(\hat{\bmbeta}-\bm{\beta})$. We derive nonasymptotic error bounds for this approximation in Section~\ref{sec:errorbounds} after the development of our method in the high-dimensional setting. If $\bm{\beta}$ is specified in the null hypothesis in a significance test, then samples from $\pi$ can be used to calculate p-values. This aspect will be explored in Section~\ref{sec:pvalue}.

\section{MCMC algorithms}\label{sec:MC}

In this section, we develop MCMC algorithms to sample from $\pi$ given $\bm{\beta}$
and $f_{\mathbf{U}}$ (or $\sigma^2$). Before that,
we first introduce a direct sampling approach which includes the residual bootstrap method as a special case.
\begin{routine}[Direct sampler]\label{rt:directsample}
Assume the error distribution is $\mathcal{D}(\mathbf{0},\sigma^2 \mathbf{I}_n)$. For $t=1,\ldots,L$
\begin{itemize}
\item[(1)]{draw $\bm{\varepsilon}^{(t)}\sim \mathcal{D}(\mathbf{0},\sigma^2 \mathbf{I}_n)$ and set $\bfy^{(t)} = \mathbf{X} \bm{\beta} + \bm{\varepsilon}^{(t)}$;}
\item[(2)]{find the minimizer $\hat{\bmbeta}^{(t)}$ of \eqref{eq:lassoloss} with $\bfy^{(t)}$ in place of $\bfy$;}
\item[(3)] if needed, calculate the subgradient vector $\mathbf{S}^{(t)}=(n\lambda\mathbf{W})^{-1}\mathbf{X}^{\mathsf{T}}(\bfy^{(t)}-\mathbf{X}\hat{\bmbeta}^{(t)})$.
\end{itemize}
\end{routine}
This approach directly draws $\bfy^{(t)}$ from its sampling distribution and requires a numerical optimization
algorithm in step (2) for each sample. Moreover, step (1) will be complicated if
we cannot draw independent samples from $\mathcal{D}(\mathbf{0},\sigma^2 \mathbf{I}_n)$. If $\bm{\varepsilon}^{(t)}$ is drawn by
resampling residuals, then Routine~\ref{rt:directsample} is equivalent to the bootstrap method of \cite{KnightFu00}.

As the density ${\pi}(\mathbf{b}_A,\mathbf{s}_I,A)$ \eqref{eq:jointdensity} has a closed-form expression given $\bm{\beta}$
and $f_{\mathbf{U}}$, MCMC and IS
can be applied to sample from and calculate expectations with respect to the distribution.
These methods may offer much more flexible and efficient alternatives to the direct sampling approach,
although the samples are either dependent or weighted.
In what follows, we propose a few special designs targeting at different applications to
exemplify the use of MCMC methods. Examples of IS will be given in Section~\ref{sec:pvalue} under
the high-dimensional setting.

\subsection{Reversibility}\label{sec:reversibility}

Our goal is to design a reversible Markov chain on the space $\Omega$, which is composed
of a finite number of subspaces $\Omega_A$, each having the same dimension $p$. Therefore, moves with an ordinary
Metropolis-Hastings (MH) ratio are sufficient, which can be seen as follows.
For any $(\mathbf{b}_A,\mathbf{s}_I,A)\in\Omega$, let $\bm{\theta}=(\theta_1,\ldots,\theta_p)$ with components given by
\begin{equation}\label{eq:definetheta}
\theta_j=\left\{\begin{array}{ll}
b_j & \mbox{ if } j\in A \\
s_j & \mbox { otherwise},
\end{array}\right.
\end{equation}
i.e., $\bm{\theta}_A=\mathbf{b}_A$ and $\bm{\theta}_I=\mathbf{s}_I$.
Then our target distribution is
$\pi(\bm{\theta}_A,\bm{\theta}_I,A)\xi_p(d\bm{\theta})$.
Suppose that $(\bm{\theta},A)$ is the current state and we
have a proposal for a new state $(\bm{\theta}^{\dag},A^{\dag})$.
In general, the proposal may only change some components of $\bm{\theta}$,
say $\theta_j$ for $j \in B\subseteq \{1,\ldots,p\}$, such that $\bm{\theta}^{\dag}_{-B}=\bm{\theta}_{-B}$. Let
$q((\bm{\theta},A),(\bm{\theta}^{\dag},A^{\dag}))$ be the density of this proposal with respect to $\xi_{|B|}$ and $I^{\dag}=\{1,\ldots,p\}\setminus A^{\dag}$. The MH ratio in terms of probability measures is
\begin{eqnarray}
& &\min\left\{1, \frac{\pi(\bm{\theta}^{\dag}_{A^{\dag}},\bm{\theta}^{\dag}_{I^{\dag}},A^{\dag})
\xi_p(d\bm{\theta}^{\dag})}{\pi(\bm{\theta}_A,\bm{\theta}_I,A)\xi_p(d\bm{\theta})} \frac{q((\bm{\theta}^{\dag},A^{\dag}),(\bm{\theta},A))
\xi_{|B|}(d\bm{\theta}_B)}{q((\bm{\theta},A),(\bm{\theta}^{\dag},A^{\dag}))\xi_{|B|}(d\bm{\theta}^{\dag}_B)} \right\} \nonumber\\
&=& \min\left\{1, \frac{\pi(\bm{\theta}^{\dag}_{A^{\dag}},\bm{\theta}^{\dag}_{I^{\dag}},A^{\dag})}{\pi(\bm{\theta}_A,\bm{\theta}_I,A)}
\frac{q((\bm{\theta}^{\dag},A^{\dag}),(\bm{\theta},A))}{q((\bm{\theta},A),(\bm{\theta}^{\dag},A^{\dag}))}
\frac{\xi_{p-|B|}(d\bm{\theta}^{\dag}_{-B})}{\xi_{p-|B|}(d\bm{\theta}_{-B})} \right\}.
\label{eq:MHratio}
\end{eqnarray}
As $\bm{\theta}^{\dag}_{-B}=\bm{\theta}_{-B}$, the dominating measures in \eqref{eq:MHratio} cancel out
and the ratio reduces to a standard MH ratio involving only densities.

Now we see that our strategy of estimator augmentation plays two roles in MCMC sampling. First,
$\mathbf{S}_\mathcal{I}$ plays the role of an auxiliary variable: The target distribution $\pi$ for $(\hat{\bmbeta}_{\mathcal{A}},\mathbf{S}_\mathcal{I},\mathcal{A})$
has a closed-form density which allows one to design an MCMC algorithm, while the distribution of interest,
that for $(\hat{\bmbeta}_{\mathcal{A}},\mathcal{A})$, is a marginal distribution of $\pi$ without a closed-form density.
Second, $\mathbf{S}_\mathcal{I}$ also plays the role
of dimension matching so that the continuous components $(\hat{\bmbeta}_{\mathcal{A}},\mathbf{S}_\mathcal{I})$ always have the
same dimension in any subspace. This eliminates the need for reversible jump MCMC \citep{Green95}. On the contrary,
if we were to sample $(\hat{\bmbeta}_{\mathcal{A}},\mathcal{A})$ (assuming a closed-form approximation to its density),
moves between two subspaces of different dimensions
would require reversible jumps, which are usually much harder to design.

\subsection{The MH Lasso sampler}\label{sec:MHLasso}

We develop an MH algorithm, called the MH Lasso sampler (MLS), with coordinate-wise update. That is,
to sequentially update each $\theta_j$, $j=1,\ldots,p$, while holding other components fixed.
Suppose the current state is $(\bm{\theta},A)$.
We design four moves to propose a new state $(\bm{\theta}^{\dag},A^{\dag})$, which are grouped into two types
according to whether $A^{\dag}=A$ or not. In the following proposals,
$\theta_j^{\dag}=b_j^{\dag}$ if $j\in A^{\dag}$ and $\theta_j^{\dag}=s_j^{\dag}$ otherwise.
\begin{definition} \label{def:proposals}
Proposals in the MLS for a given $j\in\{1,\ldots,p\}$.
\begin{itemize}
\item{Parameter-update proposals:
(P1) If $j\in A$, draw $b_j^{\dag}\sim \mathcal{N}(b_j,\tau_j^2)$.
(P2) If $j \notin A$, draw $s_j^{\dag}\sim\text{Unif}(-1,1)$. Set $A^{\dag}=A$ in both (P1) and (P2).
}
\item{Model-update proposals:
(P3) If $j\in A$, set $A^{\dag}=A\setminus \{j\}$ and draw $s_j^{\dag}\sim\text{Unif}(-1,1)$.
(P4) If $j \notin A$, set $A^{\dag}=A\cup \{j\}$ and draw $b_j^{\dag}\sim \mathcal{N}(0,\tau_j^2)$.
}
\end{itemize}
\end{definition}

The two parameter-update proposals, (P1) and (P2), are symmetric. They only change the value of $\theta_j$ and leave
$A^{\dag} =A$ so that $\det \mathbf{D}(A^{\dag})=\det \mathbf{D}(A)$. From \eqref{eq:jointdensity},
one sees that the MH ratio is simply
\begin{equation*}
\min\left\{1,\frac{f_{\mathbf{U}}(\mathbf{H}(\bm{\theta}^{\dag}_{A},\bm{\theta}^{\dag}_{I},A;\bm{\beta}))}
{f_{\mathbf{U}}(\mathbf{H}(\bm{\theta}_A,\bm{\theta}_I,A;\bm{\beta}))}\right\},
\end{equation*}
which can be computed very efficiently, especially for a normal error distribution.
The proposal (P3) removes a variable from the active set and (P4) adds a variable to the active set.
Both propose moves between two subspaces. The two proposals are the
reverse of each other and have a simple one-dimensional density.
To be concrete, the MH ratio for proposal (P3) is
\begin{equation*}
\min\left\{1, \frac{\pi(\bm{\theta}^{\dag}_{A^{\dag}},\bm{\theta}^{\dag}_{I^{\dag}},A^{\dag})}
{\pi(\bm{\theta}_A,\bm{\theta}_I,A)}\cdot\frac{\phi(b_j; 0, \tau^2_j)}{1/2} \right\},
\end{equation*}
and analogously for proposal (P4).
One needs to calculate the ratio between two determinants for these MH ratios,
\begin{equation}\label{eq:detDratio}
\frac{|\det \mathbf{D}(A^\dag)|}{|\det \mathbf{D}(A)|} = \frac{\det\mathbf{C}_{A^{\dag}A^{\dag}}}{\det\mathbf{C}_{AA}} (w_j\lambda)^{|A|-|A^{\dag}|}
\end{equation}
by \eqref{eq:defD}. As the two sets $A$ and $A^{\dag}$ differ by only one element, the ratio on the right-hand side can be calculated efficiently. When $|A|$ is large, we use the sweep operator to dynamically update $\mathbf{C}_{AA}^{-1}$ (the inverse of $\mathbf{C}_{AA}$) and obtain the ratio. See Appendix for further details. In general, however, a model-update proposal is more time-consuming than a parameter-update proposal.

This computational efficiency consideration motivates
the following scheme in the MLS which uses both types of proposals.
Let $K$ be an integer between 1 and $p$ and $\bm{\alpha}=(\alpha_j)_{1:p}$ be a vector with every
$\alpha_j>0$.
\begin{routine}[MLS]\label{rt:MHLS}
Suppose the current state is $(\bm{\theta}^{(t)},A^{(t)})$.
\begin{itemize}
\item[(1)]{Draw $K$ elements without replacement from $\{1,\ldots,p\}$ with the probability of drawing
$j$ proportional to $\alpha_j$ for each $j$. Let $M^{(t)}$ be the set of the $K$ elements.}
\item[(2)]{For $j \in M^{(t)}$, sequentially update each $\theta_j$ and the active set $A$
by an MH step with a model-update proposal.}
\item[(3)]{For $j \notin M^{(t)}$, sequentially update each $\theta_j$ by an MH step with a parameter-update proposal.}
\end{itemize}
After the above $p$ MH steps in an iteration, the state is updated to $(\bm{\theta}^{(t+1)},A^{(t+1)})$.
\end{routine}

The MLS has three input parameters, $K$, $\bm{\alpha}$, and $(\tau^2_j)_{1:p}$.
Specification of these parameters that gives good empirical performance will be provided in the numerical
examples (Section~\ref{sec:numerical}).

\subsection{The Gibbs Lasso sampler}

Let $a_j=\bm{1}(j \in A)$ and $\mathbf{a}=(a_j)_{1:p}$.
Conditional densities $\pi(\theta_j, a_j \mid \bm{\theta}_{-j}, \mathbf{a}_{-j})$ can be derived from the joint density $\pi$,
which allows for the development of a Gibbs sampler. However, as each conditional
sampling step involves calculation of one-dimensional integrals and sampling from truncated distributions,
the Gibbs sampler is more time-consuming and less efficient than the MLS for all examples on which
we have tested these algorithms.

\subsection{Conditioning on active set}

Suppose that we have constructed a Lasso-type estimate $\hat{\bmbeta}^*$ from an observed dataset
and the set of selected variables is $A^*$, which defines an estimated model.
One may want to study the sampling distribution of the estimator given the estimated model, i.e.,
$[\hat{\bmbeta}_{A^*} \mid \mathcal{A} = A^*]$. Confidence intervals of penalized estimators
have been constructed by approximating this distribution via local expansion of the $\ell_1$ norm \citep{FanLi01,Zou06}.
Since local approximation
may not be accurate for a finite sample, Monte Carlo sampling from this conditional distribution
may provide more accurate results. However, the direct sampling approach is not applicable
in practice, because $\mathcal{A} = A^*$ is often a rare event unless $p$ is very small. On the contrary, it is very efficient
to draw samples by an MH algorithm from the conditional distribution
\begin{equation}\label{eq:condAdensity}
{\pi}(\mathbf{b}_{A^*},\mathbf{s}_{I^*}\mid A^*) \propto {f}_{\mathbf{U}}(\mathbf{H}(\mathbf{b}_{A^*},\mathbf{s}_{I^*},A^*;\bm{\beta})),
\end{equation}
where $I^*=\{1,\ldots,p\}\setminus A^*$, according to \eqref{eq:condforbs}.
The distribution of interest, $[\hat{\bmbeta}_{A^*} \mid \mathcal{A} = A^*]$, is a marginal
distribution of \eqref{eq:condAdensity}.
Since evaluation of this density does not involve calculation of determinants,
each MH step is very fast.
\begin{routine}[MLS given active set]\label{rt:MLSfixedA}
Given the current state $(\mathbf{b}_{A^*}^{(t)},\mathbf{s}_{I^*}^{(t)})$,
sequentially draw $b_j^{(t+1)}$ for each $j\in A^*$ by an MH step with proposal (P1) and
$s_j^{(t+1)}$ for each $j \notin A^*$ with proposal (P2) in one iteration.
\end{routine}

\subsection{Reparameterization view}\label{sec:viewofMCMC}

To ease notation, write $\bm{\theta}=(\bm{\theta}_A,\bm{\theta}_I)$ and $\mathbf{H}_A(\bm{\theta})=\mathbf{H}(\bm{\theta},A)=\mathbf{H}(\bm{\theta}_A,\bm{\theta}_I,A;\bm{\beta})$ for $\bm{\theta}$ defined in \eqref{eq:definetheta}. Suppose that $(\bm{\theta}^{(t)},A^{(t)})$ are simulated by the MLS (Routine~\ref{rt:MHLS}) and let $\mathbf{u}^{(t)}=\mathbf{H}(\bm{\theta}^{(t)},A^{(t)})$. Since $\mathbf{H}$ is a bijection, $\mathbf{u}^{(t)}$ is a Markov chain that leaves $f_{\mathbf{U}}$ invariant. Therefore, the MLS can be understood as an MH algorithm targeting at $f_{\mathbf{U}}$ with moves designed under local reparameterization, $\bm{\theta}=\mathbf{H}_{A}^{-1}(\mathbf{u})$ for $\mathbf{u} \in U_A=\mathbf{H}_A(\Omega_A;\bm{\beta})$ \eqref{eq:subspace}. The Jacobian of this reparameterization is $[\mathbf{D}(A)]^{-1}$. Under this view, the MH ratio for a proposal $\mathbf{u}^{\dag}=\mathbf{H}(\bm{\theta}^{\dag},{A^\dag})$ given the current $\mathbf{u}$ is
\begin{equation*}
\min\left\{1, \frac{f_{\mathbf{U}}(\mathbf{u}^{\dag}) }{f_{\mathbf{U}}(\mathbf{u})}
\frac{ q((\bm{\theta}^{\dag},A^{\dag}),(\bm{\theta},A))|\det \mathbf{D}(A)|^{-1}}{q((\bm{\theta},A),(\bm{\theta}^{\dag},A^{\dag}))|\det \mathbf{D}(A^\dag)|^{-1}} \right\},
\end{equation*}
which of course coincides with \eqref{eq:MHratio}. Moreover, when $A=A^{\dag}$ such as in proposals (P1) and (P2), the Jacobian determinants cancel out as we are using the same reparameterization $\mathbf{H}_A^{-1}$ for both the proposal and the current state. Otherwise, the ratio of the Jacobian determinants accounts for the use of different reparameterizations. Clearly, if $\mathbf{u}^{(t)}\sim f_{\mathbf{U}}$ then $(\bm{\theta}^{(t)},A^{(t)})\sim \pi$.

This view provides an insight into the computational efficiency of the MLS. Under normal error assumption, $f_{\mathbf{U}}$ is the density of a multivariate normal distribution, for which a simple MH algorithm is computationally tractable and can be quite efficient. We make a comparison with the direct sampler at a conceptual level. The direct sampler draws $\mathbf{U}$ via a linear transformation of $\bm{\varepsilon}\sim\mathcal{N}_n(\mathbf{0},\sigma^2\mathbf{I}_n)$, which costs $n$ draws from a univariate normal distribution followed by a multiplication with a size $p\times n$ matrix. After that, we find $(\bm{\theta},A)=\mathbf{H}^{-1}(\mathbf{U})$ by numerical minimization due to the lack of a closed-form inverse of the mapping $\mathbf{H}$. The MLS draws $p$ ($<n$) univariate proposals in one iteration, and does not need any numerical procedure to map $\mathbf{u}^{(t)}$ back to the space $\Omega$ of $(\bm{\theta},A)$ since the moves are by design in that space already. The mapping $\mathbf{H}$ from $(\bm{\theta}^{(t)},A^{(t)})$ to $\mathbf{u}^{(t)}$ is simple and can be calculated analytically. This is fundamentally different from direct sampling which replies on a numerical procedure to find the image of each draw of $\mathbf{U}$ under the mapping $\mathbf{H}^{-1}$. The relatively time-consuming step in the MLS is calculating the ratio \eqref{eq:detDratio} when a model-update proposal is used, which can be done by at most sweeping a $|A|\times |A|$ matrix on a single position (Appendix). Owing to sparsity, $|A|$ is usually much smaller than $p$, which greatly speeds up this step. Since the target distribution in the space of $\mathbf{U}$ has a nice unimodal density, the chain $\mathbf{u}^{(t)}$ often converges fast and has low autocorrelation. Consequently, we expect to see efficiency gain over direct sampling for the same amount of computing time, which will be confirmed numerically in the next subsection.

As in the following routine, by a special initialization such that $(\bm{\theta}^{(1)},A^{(1)})\sim\pi$, the MLS can reach equilibrium in one step, which totally removes the need for burn-in iterations. This will make our method suitable for parallel computing. See Section~\ref{sec:limitation} for a more detailed discussion. However, to demonstrate the efficiency of the MLS as an independent method, we did not use Routine~\ref{rt:dirmcmc} in the numerical results.

\begin{routine}\label{rt:dirmcmc}
Draw $(\mathbf{b}^{(1)},\mathbf{s}^{(1)})$ from the direct sampler and let $(\bm{\theta}^{(1)},A^{(1)})$ be its equivalent representation. With $(\bm{\theta}^{(1)},A^{(1)})$ as the initial state, generate $(\bm{\theta}^{(t)},A^{(t)})$ for $t=2,\ldots,N$ by an MCMC algorithm targeting at $\pi$.
\end{routine}

\subsection{Numerical examples}\label{sec:numerical}

We demonstrate with numerical examples the effectiveness of the above MCMC
algorithms by comparing against the direct sampling approach. To this end,
we simulated four datasets with different combinations of $n$, $p$,
and $\sigma^2$ (Table~\ref{tab:data}). The vector of true coefficients $\bm{\beta}_0$ has 10 nonzero components,
$\beta_{0j}=1$ for $j=1,\ldots,5$ and $\beta_{0j}=-1$ for $j=6,\ldots,10$. Each row of $\mathbf{X}$ was generated independently from $\mathcal{N}_p(\mathbf{0},\bm{\Sigma}_{\mathbf{X}})$, where the diagonal and the off-diagonal elements of $\bm{\Sigma}_{\mathbf{X}}$
are 1 and $0.25$, respectively. Given the design matrix $\mathbf{X}$, the response vector $\bfy$ was drawn
from $\mathcal{N}_n(\mathbf{X} \bm{\beta}_0, \sigma^2 \mathbf{I}_n)$.

\begin{table}[ht]
\centering
\caption{Simulated datasets for MCMC\label{tab:data}}
\vspace{0.05in}
   \begin{tabular}{c|cccc}
   \hline
     Dataset	  	&  A & B & C &D  \\
     \hline
   $(n,p,\sigma^2)$ &$(500,100,1)$ & $(500,200,1)$ & $(300,100,4)$ & $(300,200,4)$ \\
   $|A^*|$ & 23 & 22 & 25 & 57 \\
  \hline
  \end{tabular}
\end{table}

The weights $w_j$ \eqref{eq:lassoloss} were set to 1 for all the following numerical results.
The Lars package by Hastie and Efron was applied to find the solution path for each dataset. The value of
$\lambda$ was chosen by minimizing the $C_p$ criterion implemented in the package,
which determined the estimated coefficients, $\hat{\bmbeta}^*=(\hat{\beta}^*_1,\ldots,\hat{\beta}^*_p)$, of a dataset.
The number of selected variables, $|A^*|$, for each dataset is given in Table~\ref{tab:data}.
We considered two types of error distributions,
the normal distribution and the elliptically symmetric distribution.
Correspondingly, we calculated $\hat{\sigma}^2$ by \eqref{eq:estsigma2} with $\check{\bmbeta}=\hat{\bmbeta}^*$
or constructed $\hat{f}_{\mathbf{U}}$ by the approach in Section~\ref{sec:asymest}.
For all the results, step (2) of the direct sampler (Routine~\ref{rt:directsample}) was implemented with the Lars package.

We first examined the performance of the MLS on sampling from the joint distribution
\eqref{eq:jointdensity} given $\hat{\bmbeta}^*$ and $\hat{\sigma}^2$ or
$\hat{f}_{\mathbf{U}}$. Let $\omega_j = \Phi(-|\hat{\beta}^*_{j}|/\zeta_j)$ for $j=1,\ldots,p$, where $\zeta_j$ is the
standard error of $\hat{\beta}^{\text{OLS}}_j$ and $\Phi$ is the cumulative distribution function of $\mathcal{N}(0,1)$.
We set $K=p/5$ and $\alpha_j \propto \omega_j + \omega_0$, where $\omega_0=\sum_j \omega_j/(5p)$
serves as a baseline weight so that each variable has a reasonable chance to be selected for model-update proposals.
See Routine~\ref{rt:MHLS} for notations. Under this setting, if the estimate $\hat{\beta}^*_{j}$
is close to zero relative to $\zeta_j$, it will have a higher chance for model-update proposals.
The $\tau_j$ used in the proposals (Definition~\ref{def:proposals})
was set to $2\zeta_j$.
The MLS was applied to each dataset 10 times independently.
Each run consisted of $L=5,500$ iterations with the first 500 as the burn-in period.
In what follows, the sampler is abbreviated as MLSn and MLSe under the normal and the
elliptically symmetric error distributions, respectively.

Figure~\ref{fig:diagnostic}(a) is the scatter plot of the samples of $\hat{\beta}_1$ and $\hat{\beta}_{50}$,
and illustrates that the distributions of some $\hat{\beta}_j$ indeed have a point mass at zero. The histogram of the subgradient $S_{50}$ is shown in Figure~\ref{fig:diagnostic}(b) with two point masses on $\pm 1$ and otherwise continuous.
Mixing of the MLS was fast, as demonstrated with two chains in Figure~\ref{fig:diagnostic}(c),
where the initial values were chosen to be about 20 standard deviations away
from each other. Figure~\ref{fig:diagnostic}(d) shows the fast decay of the autocorrelation among the samples of a $\hat{\beta}_j$, decreasing to below 0.05 in 10 to 15 iterations. The acceptance rate of the model-update proposals was generally between 0.2 and 0.4.
For the parameter-update proposals, the acceptance rate was between 0.2 and 0.4 for (P1)
and was higher than 0.6 for (P2), which is an independent proposal.

From the MCMC samples, we estimated the selection
probability $P_{s,j}=P(\hat{\beta}_j \ne 0)$, the 2.5\% and the 97.5\% quantiles
of $\hat{\beta}_j$, and the mean and the standard deviation of the conditional distribution $[\hat{\beta}_j \mid \hat{\beta}_j \ne 0]$
for each $j$. Since theoretical values are not available, we applied
the direct sampling approach to simulate 5,000 independent samples for each dataset under the normal error distribution.
These independent samples were used to estimate the above quantities as the ground truth.
The MSEs across 10 independent runs of the MLS were calculated,
and reported in Table~\ref{tab:MSEjoint} are the average MSEs over all $j$ for estimating the above five quantities.
One clearly sees that all the estimates were very accurate.
The MSE of the MLSe was greater than, but on the same order as, that of the MLSn for most estimates,
which is expected due to the loss of efficiency without assuming a normal error distribution.

\begin{table}[t]
\caption{MSE comparison for simulation from the joint sampling distribution \label{tab:MSEjoint}}
\begin{center}
   \begin{tabular}{ccccccc}
   \hline
      & Method	 & $P_{s}$	& 2.5\%	& 97.5\% & mean & SD \\
     \hline
       & MLSn & $3.38\times 10^{-4}$ &	$1.82\times 10^{-5}$	& $1.79\times 10^{-5}$ &	$4.36\times 10^{-6}$	& $2.78\times 10^{-6}$ \\
    A & MLSe & 1.29&1.20&1.19&	0.97&	1.38\\
\vspace{0.05in}  & DSn &1.11 &	2.28&	2.45& 2.23 &	2.53\\
       & MLSn & $2.13\times 10^{-4}$&	$2.89\times 10^{-5}$	&$1.74\times 10^{-5}$	& $1.22\times 10^{-5}$	& $8.44\times 10^{-6}$\\
    B & MLSe & 1.20	&	1.07 &	1.10 &	1.08 &	1.30 \\
\vspace{0.05in} & DSn & 1.26 &	1.97 &	1.89 &	2.29 &	2.74 \\
        & MLSn &  $4.14\times 10^{-4}$ &	$1.23\times 10^{-4}$&$1.24\times 10^{-4}$&	$3.20\times 10^{-5}$	& $2.28\times 10^{-5}$  \\
     C  & MLSe & 1.55 &	2.09 &	1.78 &	1.03	&2.46\\
 \vspace{0.05in}& DSn & 0.47 &	1.18 &	1.33 &	1.24 &	1.39 \\
       & MLSn & $4.34\times 10^{-4}$ & $2.96\times 10^{-4}$	& $2.85\times 10^{-4}$ &$6.37\times 10^{-5}$	& $5.02\times 10^{-5}$\\
    D & MLSe & 2.74 &		3.81 &	3.61 &	1.17 &	5.70 \\
         & DSn & 0.69 &	1.52 &	1.34 &	1.21 &	1.57 \\
     \hline
  \end{tabular}
 \end{center}
Note:  For the MLSe and the DSn, reported is the ratio of MSE to that of the MLSn. The sweep operator was used in the MLS to calculate determinant ratios for dataset D.
\end{table}

We compared the efficiency of the MLS against the direct sampler (DSn)
under the same amount of running time and under the same normal error distribution.
The DSn generated around 500 samples
in the same amount of time for 5,500 iterations of the MLSn.
The ratio of the MSE of the DSn
to that of the MLSn was calculated for each estimate (Table~\ref{tab:MSEjoint}).
For most estimates, the MLSn seems to be more efficient and may reduce the MSE
by 10\% to 60\%. The improvement was more significant for datasets A and B where the sample size $n=500$.
For the other two datasets, the MLSn showed a higher MSE in estimating selection probabilities
but was more accurate for all other estimates. Furthermore, if the error distribution is more complicated
such that one cannot simulate samples independently from the distribution, the efficiency of the direct sampler may be even lower.
These results clearly confirm the notion that the MLS can serve as an efficient alternative to the direct sampling method
for simulating from the sampling distribution of a Lasso-type estimator.

Next, we implemented Routine~\ref{rt:MLSfixedA} to sample from the conditional distribution of $\hat{\bmbeta}$ given the model
selected according to the $C_p$ criterion, i.e., $[\hat{\bmbeta}_{A^*} \mid \mathcal{A} =A^*]$ with $|A^*|$ given in Table~\ref{tab:data}.
The same parameter setting as that in the previous example was used to run the MLSn and the MLSe.
We estimated the 2.5\% and the 97.5\% quantiles, the mean, and the standard deviation of $\hat{\beta}_j$
for $j \in A^*$. The model space is composed of
$2^p$ models, and the probability of the model $A^*$, $P(\mathcal{A} =A^*)$, is practically zero
for the datasets used here. Therefore, the direct sampling approach is not applicable.
This shows the advantage and flexibility of the Monte Carlo algorithms.
Since we cannot construct ground truth for this example, the accuracy of an estimate is measured by its variance across 10 independent runs of the MLS, averaging over $j\in A^*$ (Table~\ref{tab:MSEcond}).
The variance of every estimate was on the order of $10^{-5}$ or smaller for datasets A, B and C and was
on the order of $10^{-4}$ or smaller for dataset D under both error models. This highlights the stability of
the MLS in approximating sampling distributions across different runs. There were cases in which
the variance of the MLSe was smaller. This does not necessarily suggest that the MLSe provided a more accurate
estimate, as the loss of efficiency without the normal error assumption is likely to result in a higher bias.

\begin{table}[t]
\centering
\caption{Variance comparison for conditional sampling given active set\label{tab:MSEcond}}
\vspace{0.05in}
   \begin{tabular}{cccccc}
   \hline
      & Method	 & 2.5\%	& 97.5\% & mean & SD \\
     \hline
   				A    & MLSn & $1.21\times 10^{-5}$& $1.28\times 10^{-5}$	& $2.21\times 10^{-6}$	& $1.03\times 10^{-6}$ \\
 \vspace{0.05in}    & MLSe &0.90 &	1.02 &	0.92 &	1.05 \\

   				B    & MLSn & $1.47\times 10^{-5}$ &	$1.19\times 10^{-5}$	& $3.19\times 10^{-6}$ &	$9.60\times 10^{-7}$ \\
 \vspace{0.05in}    & MLSe &1.22 &	1.15 &	1.02 &	1.23 \\

 				C   & MLSn & $7.66\times 10^{-5}$	 & $8.65\times 10^{-5}$	 &$1.59\times 10^{-5}$& $7.08\times 10^{-6}$ \\
 \vspace{0.05in}    & MLSe &1.07 &	0.95 &	1.01 &	1.00 \\

 				D & MLSn & $1.67\times 10^{-4}$	&$1.78\times 10^{-4}$ &	$2.55\times 10^{-5}$	 &$1.28\times 10^{-5}$\\
   & MLSe &	0.77 &	0.97 &	0.77 &	0.79 \\
     \hline
  \end{tabular}

\vspace{0.05in}Note: Variance of the MLSe is reported as the ratio to that of the MLSn.
\end{table}

\section{High-dimensional setting}\label{sec:highdim}

Recent efforts have established theoretical properties of $\ell_1$-penalized linear regression
in high dimension with $p>n$ \citep{Meinshausen06,ZhaoYu06,ZhangHuang08,Bickel09}.
Under this setting, we assume $\text{rank}(\mathbf{X})=n<p$. Consequently, $\text{row}(\mathbf{X})$ is an $n$-dimensional subspace of $\mathbb{R}^p$ and the Gram matrix $\mathbf{C}$ has $n$ positive eigenvalues, denoted by $\Lambda_j>0$, $j=1,\ldots,n$. The associated orthonormal eigenvectors $\mathbf{v}_j \in \mathbb{R}^p$, $j=1,\ldots,n$, form a basis for $\text{row}(\mathbf{X})$.
Choose orthonormal vectors $\mathbf{v}_{p+1},\ldots,\mathbf{v}_p$ to form a basis for the null space of $\mathbf{X}$, $\text{null}(\mathbf{X})$, and let $\mathbf{V}=(\mathbf{v}_1 | \ldots |\mathbf{v}_p)$. Then $R=\{1,\ldots,n\}$ and $N=\{n+1,\ldots,p\}$ index the columns of $\mathbf{V}$ that form respective bases for $\text{row}(\mathbf{X})$ and $\text{null}(\mathbf{X})$.

\begin{assumption}\label{as:X}
Every $n$ columns of $\mathbf{X}$ are linearly independent and every $(p-n)$ rows of $\mathbf{V}_N$ are linearly independent.
\end{assumption}
The first part of this assumption is sufficient for the columns of $\mathbf{X}$ being in general position, which guarantees that $(\hat{\bmbeta},\mathbf{S})$ is unique for any $\bfy$ and $\lambda>0$ (Lemma~\ref{lm:prelim}). The second part will ease our derivation of the joint density of the augmented estimator. Note that Assumption~\ref{as:X} holds with probability one if the entries of $\mathbf{X}$ are drawn from a continuous distribution on $\mathbb{R}^{n\times p}$.

\subsection{The bijection}

Although $\mathbf{U}=\frac{1}{n}\mathbf{X}^{\mathsf{T}}\bm{\varepsilon}$ is a $p$-vector, by definition it always lies in $\text{row}(\mathbf{X})$.
Therefore, $\mathbf{V}_N^{\mathsf{T}} \mathbf{U}=\mathbf{0}$ and the $n$-vector $\mathbf{R}=\mathbf{V}_R^{\mathsf{T}}\mathbf{U}$
gives the coordinates of $\mathbf{U}$ with respect to the basis $\mathbf{V}_R$.
If $\bm{\varepsilon}$ follows a continuous distribution on $\mathbb{R}^n$, then $\mathbf{R}$ has
a proper density with respect to $\xi_n$. For example, if $\bm{\varepsilon}\sim \mathcal{N}_n(\mathbf{0},\sigma^2 \mathbf{I}_n)$,
then $\mathbf{R}\sim\mathcal{N}_n(\mathbf{0},\frac{\sigma^2}{n} \bm{\Lambda})$ with $\bm{\Lambda}=\mbox{diag}(\Lambda_1,\ldots,\Lambda_n)$. Now $\mathbf{R}$ plays the same role
as $\mathbf{U}$ does in the low-dimensional case. We will use the known distribution of $\mathbf{R}$
to derive the distribution of the augmented estimator $(\hat{\bmbeta}_\mathcal{A},\mathbf{S}_\mathcal{I},\mathcal{A})$.

However, a technical difficulty is that when $p>n$, the map $\mathbf{H}$ defined in \eqref{eq:lassojoint}
is not a mapping from $\Omega$
to $\text{row}(\mathbf{X})$ as $\mathbf{H}(\mathbf{b}_A,\mathbf{s}_I,A;\bm{\beta})\in\mathbb{R}^p$ is not necessarily in $\text{row}(\mathbf{X})$
for every $(\mathbf{b}_A,\mathbf{s}_I,A)\in \Omega$. We thus need to remove those ``illegal"
points in $\Omega$ so that the image of $\mathbf{H}$ always lies in the row space of $\mathbf{X}$. This is achieved by
imposing the constraint that
\begin{equation}\label{eq:constraintinU}
\mathbf{V}_N^{\mathsf{T}} \mathbf{H}(\hat{\bmbeta}_{\mathcal{A}},\mathbf{S}_\mathcal{I},\mathcal{A};\bm{\beta})=\mathbf{V}_N^{\mathsf{T}} \mathbf{U}=\mathbf{0},
\end{equation}
i.e., the image of $\mathbf{H}$ must be orthogonal to $\text{null}(\mathbf{X})$. It is more convenient to use the equivalent definition of $\mathbf{H}$
in \eqref{eq:lassoKKTinU}, i.e.,
\begin{equation}\label{eq:equivH}
\mathbf{H}(\hat{\bmbeta}_{\mathcal{A}},\mathbf{S}_\mathcal{I},\mathcal{A};\bm{\beta})=\mathbf{C} \hat{\bmbeta} + \lambda \mathbf{W} \mathbf{S} - \mathbf{C} \bm{\beta}.
\end{equation}
Because $\mathbf{C}(\hat{\bmbeta}-\bm{\beta})=\frac{1}{n}\mathbf{X}^{\mathsf{T}}\mathbf{X}(\hat{\bmbeta}-\bm{\beta})\in\text{row}(\mathbf{X})$,
constraint \eqref{eq:constraintinU} is equivalent to
\begin{equation}\label{eq:constraintonS}
 \mathbf{V}_N^{\mathsf{T}} \mathbf{W} \mathbf{S} = \mathbf{V}_{\mathcal{A} N}^{\mathsf{T}} \mathbf{W}_{\mathcal{A}\calA} \mathbf{S}_\mathcal{A}
 + \mathbf{V}_{\mathcal{I} N}^{\mathsf{T}} \mathbf{W}_{\mathcal{I}\calI} \mathbf{S}_\mathcal{I} = \mathbf{0}.
\end{equation}
In words, the constraint is that
the vector $\mathbf{W}\mathbf{S}$ must lie in $\text{row}(\mathbf{X})$. Therefore, we have a more
restricted space for the augmented estimator $(\hat{\bmbeta}_{\mathcal{A}},\mathbf{S}_\mathcal{I},\mathcal{A})$ in the high-dimensional case,
\begin{equation}\label{eq:spacehigh}
\Omega_r=\{(\mathbf{b}_A,\mathbf{s}_I,A)\in \Omega: \mathbf{V}_{AN}^{\mathsf{T}} \mathbf{W}_{AA} \mbox{sgn}(\mathbf{b}_A)
 + \mathbf{V}_{IN}^{\mathsf{T}} \mathbf{W}_{II} \mathbf{s}_I =\mathbf{0} \}.
\end{equation}
Restricted to this space, $\mathbf{H}$ is a bijection.
\begin{lemma}\label{lm:bijectioninhighdim}
If $p>n$ and Assumption~\ref{as:X} holds, then for any $\bm{\beta}$ and $\lambda>0$, the restriction of the mapping $\mathbf{H}$ \eqref{eq:lassojoint} to $\Omega_r$, denoted by $\mathbf{H}\mid_{\Omega_r}$,
is a bijection that maps $\Omega_r$ onto $\text{row}(\mathbf{X})$.
\end{lemma}
\begin{proof}
For any $\mathbf{U} \in \text{row}(\mathbf{X})$, there is a unique $(\hat{\bmbeta}_\mathcal{A},\mathbf{S}_\mathcal{I},\mathcal{A})$
such that $\mathbf{U}=\mathbf{H}(\hat{\bmbeta}_{\mathcal{A}},\mathbf{S}_\mathcal{I},\mathcal{A};\bm{\beta})\in\text{row}(\mathbf{X})$ by Lemma~\ref{lm:prelim}.
Thus, $(\hat{\bmbeta}_\mathcal{A},\mathbf{S}_\mathcal{I},\mathcal{A})$ satisfies the constraint
\eqref{eq:constraintinU} and lies in $\Omega_r$.
For any $(\hat{\bmbeta}_\mathcal{A},\mathbf{S}_\mathcal{I},\mathcal{A}) \in \Omega_r$,
$\mathbf{V}_N^{\mathsf{T}} \mathbf{H}(\hat{\bmbeta}_{\mathcal{A}},\mathbf{S}_\mathcal{I},\mathcal{A};\bm{\beta})=\mathbf{0}$
and $\mathbf{H}$ maps it into $\text{row}(\mathbf{X})$.
\end{proof}
\begin{remark}\label{rm:dimofH}
Fixing $\mathcal{A}=A$, \eqref{eq:constraintonS} specifies $|N|=p-n$ constraints,
and thus, the continuous components  $(\hat{\bmbeta}_A,\mathbf{S}_I)\in\mathbb{R}^p$
lie in an $n$-dimensional subspace of $\Omega_A$ \eqref{eq:subspace}.
The bijection $\mathbf{H}\mid_{\Omega_r}$
maps a finite number of $n$-dimensional subspaces onto $\text{row}(\mathbf{X})$ which is an $\mathbb{R}^n$.
\end{remark}
Now we represent the bijection $\mathbf{H}\mid_{\Omega_r}$ in terms of its coordinates with respect to $\mathbf{V}_R$
and equate it with $\mathbf{R}=\mathbf{V}_R^{\mathsf{T}}\mathbf{U}$:
\begin{equation}\label{eq:HinR}
\mathbf{R}=\mathbf{V}_R^{\mathsf{T}} \mathbf{H}(\hat{\bmbeta}_\mathcal{A},\mathbf{S}_\mathcal{I},\mathcal{A};\bm{\beta}) \defi \mathbf{H}_r(\hat{\bmbeta}_\mathcal{A},\mathbf{S}_\mathcal{I},\mathcal{A};\bm{\beta}).
\end{equation}

\subsection{Joint sampling distribution}

The distribution for $(\hat{\bmbeta}_{\mathcal{A}},\mathbf{S}_{\mathcal{I}},\mathcal{A}) \in \Omega_r$
is completely given by the distribution of $\mathbf{R}$ via the bijective map $\mathbf{H}_r:\Omega_r \to \mathbb{R}^n$.
The only task left is to determine the Jacobian of $\mathbf{H}_r$, taking into account the constraint \eqref{eq:constraintonS}.
Left Multiplying by $\mathbf{V}_R^{\mathsf{T}}$ both sides of Equation~\eqref{eq:lassojoint_a},
with the simple facts that $\mathbf{V}_R^{\mathsf{T}}\mathbf{W}_{\mathcal{A}}=\mathbf{V}_{\mathcal{A} R}^{\mathsf{T}}\mathbf{W}_{\mathcal{A}\calA}$ and
$\mathbf{V}_R^{\mathsf{T}}\mathbf{W}_{\mathcal{I}}=\mathbf{V}_{\mathcal{I} R}^{\mathsf{T}} \mathbf{W}_{\mathcal{I}\calI}$, gives
\begin{equation}\label{eq:Hinhigh}
\mathbf{R} =  \mathbf{V}_R^{\mathsf{T}} \mathbf{C}_{\mathcal{A}} \hat{\bmbeta}_{\mathcal{A}} + \lambda \mathbf{V}_{\mathcal{A} R}^{\mathsf{T}}\mathbf{W}_{\mathcal{A}\calA} \mathbf{S}_{\mathcal{A}}
+\lambda \mathbf{V}_{\mathcal{I} R}^{\mathsf{T}} \mathbf{W}_{\mathcal{I}\calI} \mathbf{S}_{\mathcal{I}} -\mathbf{V}_R^{\mathsf{T}}\mathbf{C}\bm{\beta}.
\end{equation}
For any fixed value of $\mathcal{A}$, differentiating $\mathbf{R}$ and both sides of the constraint \eqref{eq:constraintonS}
with respect to $(\hat{\bmbeta}_{\mathcal{A}},\mathbf{S}_\mathcal{I})$ give, respectively,
\begin{eqnarray}
& d\mathbf{R} =  \mathbf{V}_R^{\mathsf{T}} \mathbf{C}_{\mathcal{A}} d \hat{\bmbeta}_{\mathcal{A}}
+\lambda \mathbf{V}_{\mathcal{I} R}^{\mathsf{T}} \mathbf{W}_{\mathcal{I}\calI} d \mathbf{S}_{\mathcal{I}}, \label{eq:dR}\\
&\mathbf{V}_{\mathcal{I} N}^{\mathsf{T}} \mathbf{W}_{\mathcal{I}\calI} d \mathbf{S}_{\mathcal{I}}=\mathbf{0}. \label{eq:dS}
\end{eqnarray}
Therefore, the constraint implies that $d \mathbf{S}_{\mathcal{I}}$ is in $\text{null} (\mathbf{V}_{\mathcal{I} N}^{\mathsf{T}} \mathbf{W}_{\mathcal{I}\calI})$.
\begin{lemma}\label{lm:rank}
If $p>n$ and Assumption~\ref{as:X} is satisfied, then the dimension of $\text{null} (\mathbf{V}_{\mathcal{I} N}^{\mathsf{T}} \mathbf{W}_{\mathcal{I}\calI})$ is $n-|\mathcal{A}|\geq 0$.
\end{lemma}
\begin{proof}
Under the assumption, the minimizer $\hat{\bmbeta}$ of \eqref{eq:lassoloss} is unique and always has an active set with size
$|\mathcal{A}|\leq \min\{n,p\}=n$. See Lemma 14 in \cite{Tibshirani13}.
If Assumption~\ref{as:X} is satisfied, any $|N|$ rows of $\mathbf{V}_N$ are linearly independent. Since $|\mathcal{I}|=p-|\mathcal{A}|\geq p-n=|N|$,
the rank of the $|N| \times |\mathcal{I}|$ matrix, $\mathbf{V}_{\mathcal{I} N}^{\mathsf{T}} \mathbf{W}_{\mathcal{I}\calI}$, is $p-n$.
Then it follows that the dimension of $\text{null} (\mathbf{V}_{\mathcal{I} N}^{\mathsf{T}} \mathbf{W}_{\mathcal{I}\calI})$ is $|\mathcal{I}|-(p-n)=n-|\mathcal{A}|\geq 0$.
\end{proof}
Let $\mathbf{B}(\mathcal{I})\in\mathbb{R}^{|\mathcal{I}|\times(n-|\mathcal{A}|)}$ be an orthonormal basis for $\text{null} (\mathbf{V}_{\mathcal{I} N}^{\mathsf{T}} \mathbf{W}_{\mathcal{I}\calI})$.
Let $d\tilde{\mathbf{S}}$ be the coordinates of $d\mathbf{S}_{\mathcal{I}}$ with respect to the basis $\mathbf{B}(\mathcal{I})$, i.e.,
$d\mathbf{S}_{\mathcal{I}}=\mathbf{B}(\mathcal{I})d\tilde{\mathbf{S}}$. Note that $d\tilde{\mathbf{S}} \in \mathbb{R}^{n-|\mathcal{A}|}$ according to
the above lemma. Then \eqref{eq:dR} becomes
\begin{eqnarray}
d\mathbf{R} & =  & \mathbf{V}_R^{\mathsf{T}} \mathbf{C}_{\mathcal{A}} d \hat{\bmbeta}_{\mathcal{A}}
+\lambda \mathbf{V}_{\mathcal{I} R}^{\mathsf{T}} \mathbf{W}_{\mathcal{I}\calI} \mathbf{B}(\mathcal{I}) d\tilde{\mathbf{S}}  \nonumber \\
& = & \mathbf{T}(\mathcal{A}) \left(\begin{array}{c} d \hat{\bmbeta}_{\mathcal{A}} \\ d\tilde{\mathbf{S}} \end{array}\right), \label{eq:Jacob}
\end{eqnarray}
where $\mathbf{T}(\mathcal{A})=(\mathbf{V}_R^{\mathsf{T}} \mathbf{C}_{\mathcal{A}}\mid \lambda \mathbf{V}_{\mathcal{I} R}^{\mathsf{T}} \mathbf{W}_{\mathcal{I}\calI} \mathbf{B}(\mathcal{I}))$,
an $n\times n$ matrix, is the Jacobian of the map $\mathbf{H}_r$.
The dimension of $(d \hat{\bmbeta}_{\mathcal{A}},d\tilde{\mathbf{S}})$ is always $n$ for any $\mathcal{A}$.
This confirms the notion in Remark~\ref{rm:dimofH} that the continuous components
$(\hat{\bmbeta}_{\mathcal{A}},\mathbf{S}_{\mathcal{I}})$ lie in an $n$-dimensional subspace when $\mathcal{A}$ is fixed.

Now we are ready to
derive the density for $(\hat{\bmbeta}_\mathcal{A},\mathbf{S}_{\mathcal{I}},\mathcal{A})$ in high dimension. For $(\mathbf{b}_A,\mathbf{s}_I,A)\in\Omega_r$,
$d\mathbf{s}_I=\mathbf{B}(I)d\tilde{\mathbf{s}}$ for some $d\tilde{\mathbf{s}} \in \mathbb{R}^{n-|A|}$ and $\xi_{n-|A|}(d\tilde{\mathbf{s}})$
gives the infinitesimal volume at $\mathbf{s}_I$ subject to constraint \eqref{eq:spacehigh}.

\begin{theorem}\label{thm:jointhigh}
Assume that $p>n$ and Assumption~\ref{as:X} holds. Let $f_{\mathbf{R}}$ be the probability density
of $\mathbf{R}$ with respect to $\xi_n$. For $(\mathbf{b}_A,\mathbf{s}_I,A)\in\Omega_r$,
the joint distribution of $(\hat{\bmbeta}_\mathcal{A},\mathbf{S}_\mathcal{I}, \mathcal{A})$ is given by
\begin{eqnarray}\label{eq:jointdensityhigh}
 P(\hat{\bmbeta}_A\in d\mathbf{b}_A,\mathbf{S}_I \in d\mathbf{s}_I,\mathcal{A}=A)
&=&f_{\mathbf{R}}(\mathbf{H}_r(\mathbf{b}_A,\mathbf{s}_I,A;\bm{\beta}))|\det\mathbf{T}(A)| \xi_n(d\mathbf{b}_A d\tilde{\mathbf{s}}) \nonumber \\
& \defi & \pi_r(\mathbf{b}_A,\mathbf{s}_I,A) \xi_n(d\mathbf{b}_A d\tilde{\mathbf{s}}).
\end{eqnarray}
Particularly, if $\bm{\varepsilon}\sim \mathcal{N}_n(\mathbf{0},\sigma^2 \mathbf{I}_n)$, then
\begin{equation}\label{eq:densitynormalhigh}
\pi_r(\mathbf{b}_A,\mathbf{s}_I,A)=\phi_n\left(\mathbf{H}_r(\mathbf{b}_A,\mathbf{s}_I,A;\bm{\beta});\mathbf{0},n^{-1}{\sigma^2} \bm{\Lambda}\right)|\det\mathbf{T}(A)|.
\end{equation}
\end{theorem}
\begin{proof}
The proof is analogous to that of Theorem~\ref{thm:joint}.
Let $\mathbf{r}=\mathbf{H}_r(\mathbf{b}_A,\mathbf{s}_I,A;\bm{\beta}) \in \mathbb{R}^n$. For any fixed $A$,
\begin{equation*}
d\mathbf{r} = \mathbf{T}(A) \left(\begin{array}{c} d \mathbf{b}_A \\ d\tilde{\mathbf{s}} \end{array}\right)
\end{equation*}
from \eqref{eq:Jacob} and thus $\xi_n(d\mathbf{r})=|\det\mathbf{T}(A)| \xi_n(d\mathbf{b}_A d\tilde{\mathbf{s}})$.
With the bijective nature of $\mathbf{H}_r$ and its restriction to any $A$, a change of variable gives
\begin{align*}
& P(\hat{\bmbeta}_A\in d\mathbf{b}_A,\mathbf{S}_I \in d\mathbf{s}_I,\mathcal{A}=A) =  P(\mathbf{R} \in d\mathbf{r}) \\
&\quad\quad = f_{\mathbf{R}}(\mathbf{H}_r(\mathbf{b}_A,\mathbf{s}_I,A;\bm{\beta}))|\det\mathbf{T}(A)| \xi_n(d\mathbf{b}_A d\tilde{\mathbf{s}}).
\end{align*}
If $\bm{\varepsilon}\sim \mathcal{N}_n(\mathbf{0},\sigma^2 \mathbf{I}_n)$
then $\mathbf{R}\sim\mathcal{N}_n(\mathbf{0},\frac{\sigma^2}{n} \bm{\Lambda})$, which leads to \eqref{eq:densitynormalhigh} immediately.
\end{proof}
\begin{remark}\label{rm:welldefine}
The density $\pi_r(\mathbf{b}_A,\mathbf{s}_I,A)$ does not depend on which orthonormal basis we choose
for $\text{null}(\mathbf{V}_{I N}^{\mathsf{T}} \mathbf{W}_{II})$. If $\mathbf{B}'(I)$ is another orthonormal basis, then $\mathbf{B}'(I)=\mathbf{B}(I)\mathbf{O}$,
where $\mathbf{O}$ is an $(n-|A|) \times (n-|A|)$ orthogonal matrix and $|\det\mathbf{O}|=1$.
Correspondingly,
\begin{equation*}
\mathbf{T}'(A)=(\mathbf{V}_R^{\mathsf{T}} \mathbf{C}_{A}\mid \lambda \mathbf{V}_{I R}^{\mathsf{T}} \mathbf{W}_{II} \mathbf{B}'(I))=\mathbf{T}(A)\mbox{diag}(\mathbf{I}_{|A|},\mathbf{O}),
\end{equation*}
and thus $|\det \mathbf{T}'(A)|=|\det \mathbf{T}(A)|$.
\end{remark}
\begin{remark}\label{rm:unified}
One may unify Theorems~\ref{thm:joint} and \ref{thm:jointhigh} with the use of cumbersome notations, but the idea is simple.
Note that $\mathbf{T}(A)$ and $\mathbf{D}(A)$ in \eqref{eq:lassojoint} are connected by
\begin{equation*}
\mathbf{T}(A)=\mathbf{V}_R^{\mathsf{T}} (\mathbf{C}_{A}\mid \lambda \mathbf{W}_{I}\mathbf{B}(I))=\mathbf{V}_R^{\mathsf{T}}\mathbf{D}(A) \mbox{diag}(\mathbf{I}_{|A|},\mathbf{B}(I)).
\end{equation*}
If $\text{rank}(\mathbf{X})=p\leq n$, the set $N$ reduces to the empty set and $\mathbf{V}_R=\mathbf{V}$.
Hence, the constraint \eqref{eq:constraintonS} no long exists,
the space $\Omega_r$ is the same as $\Omega$,
and $\text{null}(\mathbf{V}_{I N}^{\mathsf{T}} \mathbf{W}_{II})$ is simply $\mathbb{R}^{|I|}$ for any $I=\{1,\ldots,p\}\setminus A$.
Choosing $\mathbf{B}(I)=\mathbf{I}_{|I|}$ leads to $d\mathbf{s}_I=d\tilde{\mathbf{s}}$,
which shows that the probability in \eqref{eq:jointdensityhigh} reduces to
that in \eqref{eq:jointdensity}. In this case,
$\mathbf{T}(A)=\mathbf{V}^{\mathsf{T}}\mathbf{D}(A)$,
i.e., a column of $\mathbf{T}(A)$ gives the coordinates of the corresponding column of $\mathbf{D}(A)$
with respect to the basis $\mathbf{V}$.
\end{remark}

In principle, one can develop MCMC algorithms to sample from the joint distribution \eqref{eq:jointdensityhigh}.
Development of such an algorithm is a little tedious because of the constraints in the definition of $\Omega_r$
and the use of different bases $\mathbf{B}(I)$ in different subspaces.
However, the explicit density given in Theorem~\ref{thm:jointhigh} allows us to
develop very efficient IS algorithms for approximating tail probabilities with respect
to the sampling distribution of a Lasso-type estimator.

\section{P-value calculation by IS}\label{sec:pvalue}

To simplify description, we focus on the high-dimensional setting
with normal errors so that the distribution of interest is $\pi_r$ \eqref{eq:densitynormalhigh}.
For a fixed $\mathbf{X}$, the density $\pi_r$, the bijection $\mathbf{H}_r$ \eqref{eq:HinR} and the matrix $\mathbf{T}$ \eqref{eq:Jacob} are written
as $\pi_r(\mathbf{b}_A,\mathbf{s}_I,A;\bm{\beta},\sigma^2,\lambda)$, $\mathbf{H}_r(\mathbf{b}_A,\mathbf{s}_I,A; \bm{\beta}, \lambda)$,
and $\mathbf{T}(A; \lambda)$, respectively, to explicitly indicate their dependency on different parameters.
Suppose we are given a Lasso-type estimate $\hat{\bmbeta}^{*}$ for an
observed dataset with a tuning parameter $\lambda^{*}$. Under the null model
$\mathcal{H}: \bm{\beta}=\bm{\beta}_0, \sigma^2=\sigma^2_0$, we want to
calculate the p-value of some test statistic $T(\hat{\bmbeta}) \in \mathbb{R}$ constructed from
the Lasso-type estimator $\hat{\bmbeta}$ for $\lambda=\lambda^*$. Precisely, the desired p-value is
\begin{equation}\label{eq:pvalue}
q^*=P(|T(\hat{\bmbeta})|\geq T^*; \mathcal{H}, \lambda^*)=\int_{\Omega_r^*} \pi_r(\mathbf{b}_A,\mathbf{s}_I,A;\bm{\beta}_0,\sigma_0^2,\lambda^*)
\xi_n(d\mathbf{b}_A d\tilde{\mathbf{s}}),
\end{equation}
where $T^*=|T(\hat{\bmbeta}^*)|$ and $\Omega_r^*=\{(\mathbf{b}_A,\mathbf{s}_I,A)\in \Omega_r: |T(\mathbf{b})|\geq T^*\}$.
Even if we can directly sample from $\pi_r(\bullet;\bm{\beta}_0,\sigma_0^2,\lambda^*)$,
estimating $q^*$ will be extremely difficult when it is very small. With the closed-form density $\pi_r$,
we can use IS to solve this challenging problem.

\subsection{Importance sampling}

Our target distribution is $\pi_r(\bullet;\bm{\beta}_0,\sigma_0^2,\lambda^*)$ and
we propose to use $\pi_r(\bullet;\bm{\beta}_0,(\sigma^2)^{\dag},\lambda^{\dag})$ as a trial distribution
to estimate expectations with respect to the target distribution via IS.
First, note that the trial and the target distributions have the same support as the constraint in \eqref{eq:spacehigh}
that defines the space $\Omega_r$ only depends on $\mathbf{X}$. Thus, a sample from
the trial distribution $\pi_r(\bullet;\bm{\beta}_0,(\sigma^2)^{\dag},\lambda^{\dag})$ also satisfies the
constraint for the target distribution. Second, one can easily simulate from the trail distribution
by the direct sampler (Routine~\ref{rt:directsample}).
Third, the importance weight for a sample $(\mathbf{b}_A,\mathbf{s}_I,A)$ from the trial distribution can be calculated efficiently.
Let $(r_1(\bm{\beta},\lambda),\cdots,r_n(\bm{\beta},\lambda))=\mathbf{H}_r(\mathbf{b}_A,\mathbf{s}_I,A;\bm{\beta},\lambda)\in \mathbb{R}^n$ and note that $\det\mathbf{T}(A;\lambda)=\lambda^{n-|A|} \det \mathbf{T}(A;1)$.
Using the fact that $\bm{\Lambda}=\mbox{diag}(\Lambda_1,\ldots,\Lambda_n)$, the importance weight
\begin{eqnarray}\label{eq:ISweight}
& &\frac{\phi_n\left(\mathbf{H}_r(\mathbf{b}_A,\mathbf{s}_I,A;\bm{\beta}_0, \lambda^*);\mathbf{0},n^{-1}{\sigma_0^2} \bm{\Lambda}\right)|\det\mathbf{T}(A;\lambda^*)|}
{\phi_n\left(\mathbf{H}_r(\mathbf{b}_A,\mathbf{s}_I,A;\bm{\beta}_0, \lambda^{\dag});\mathbf{0},n^{-1}{(\sigma^2)^{\dag}} \bm{\Lambda}\right)|\det\mathbf{T}(A;\lambda^{\dag})|} \nonumber\\
&\propto &
\exp\left[\frac{n}{2(\sigma^2)^{\dag}}\sum_{i=1}^n\frac{r_i^2(\bm{\beta}_0,\lambda^{\dag})}{\Lambda_i}
-\frac{n}{2 \sigma_0^2}\sum_{i=1}^n\frac{r_i^2(\bm{\beta}_0,\lambda^{*})}{\Lambda_i} \right]
\left(\frac{\lambda^*}{\lambda^{\dag}}\right)^{n-|A|}  \nonumber \\
&\defi & w(\mathbf{b}_A,\mathbf{s}_I,A;\sigma_0^2,\lambda^*).
\end{eqnarray}
Essentially, for each sample, we only need to compute the image of the map $\mathbf{H}_r$ and two sums of squares.
\begin{routine}\label{rt:IS}
Draw $(\mathbf{b}_A,\mathbf{s}_I,A)^{(t)}$, $t=1,\ldots,L$, from the trial distribution
$\pi_r(\bullet;\bm{\beta}_0,(\sigma^2)^{\dag},\lambda^{\dag})$ by Routine~\ref{rt:directsample}.
Then the IS estimate for the p-value $q^*$ is given by
\begin{equation}\label{eq:ISest}
\hat{q}^{\text{(IS)}}=\frac{\sum_{t=1}^L w((\mathbf{b}_A,\mathbf{s}_I,A)^{(t)};\sigma_0^2,\lambda^*)\bm{1}(|T(\mathbf{b}^{(t)})|\geq T^*)}
{\sum_{t=1}^L w((\mathbf{b}_A,\mathbf{s}_I,A)^{(t)};\sigma_0^2,\lambda^*)}.
\end{equation}
\end{routine}

The key is to choose the parameters $(\sigma^2)^{\dag}$ and $\lambda^{\dag}$ in the trial distribution so that we have
a substantial fraction of samples for which $|T(\mathbf{b}^{(t)})|\geq T^*$. Next we discuss some guidance on
tuning these parameters.

\subsection{Tuning trial distributions}

We illustrate our procedure for tuning the trial distribution assuming $\bm{\beta}_0=\mathbf{0}$, i.e., the null hypothesis is
$\mathcal{H}_0: \bm{\beta}=\mathbf{0}, \sigma^2=\sigma^2_0$. In this case, the problem is difficult when
$P(\hat{\bmbeta}=\mathbf{0})$ is close to one under the target distribution $\pi_r(\bullet; \mathbf{0}, \sigma_0^2, \lambda^*)$.
In other words, $\lambda^*$ is too big to obtain any nonzero estimate of the coefficients and consequently
the p-value $P(|T(\hat{\bmbeta})|\geq T^*; \mathcal{H}_0, \lambda^*)$ becomes a tail probability.
Thus, one may want to choose the trial distribution
$\pi_r(\bullet; \mathbf{0}, (\sigma^2)^{\dag}, \lambda^{\dag})$ under which there is a higher probability for
nonzero $\hat{\bmbeta}$. In general, we achieve this by choosing $(\sigma^2)^{\dag}=M^{\dag} \sigma_0^2$ ($M^{\dag}>1$)
and then tuning $\lambda^{\dag}$ accordingly. When we increase $\sigma^2$, the variance of $\mathbf{U}$
increases and thus $\mathbf{U}$ will have a wider spread in $\text{row}(\mathbf{X})$. This will increase
the variance of the augmented estimator $(\hat{\bmbeta}_\mathcal{A},\mathbf{S}_\mathcal{I},\mathcal{A})$.
As illustrated in Figure~\ref{fig:bijection}, a larger variance in $\mathbf{U}$ will lead to a more uniform distribution over
different subspaces $\{\Omega_A\}$.

The following simple procedure is used to determine $\lambda^{\dag}$ given $(\sigma^2)^{\dag}$, which
works very well based on our empirical study.
\begin{routine}\label{rt:tuning}
Draw $\bfy^{(t)}$ from $\mathcal{N}_n(\mathbf{0},(\sigma^2)^{\dag} \mathbf{I}_n)$
and calculate
$\lambda^{(t)}=n^{-1}\|\mathbf{W}^{-1}\mathbf{X}^{\mathsf{T}}\bfy^{(t)}\|_{\infty}$
for $t=1,\ldots,L_{\text{pilot}}$. Then set $\lambda^{\dag}$ to the first quartile of
$\{\lambda^{(t)}: t=1,\ldots,L_{\text{pilot}}\}$.
\end{routine}
Setting $\hat{\bmbeta}=\mathbf{0}$ in \eqref{eq:lassograd}, we have
\begin{equation*}
n^{-1}\|\mathbf{W}^{-1}\mathbf{X}^{\mathsf{T}}\bfy\|_{\infty}=\lambda \|\mathbf{S}\|_{\infty}\leq \lambda,
\end{equation*}
which shows that the $\lambda^{(t)}$ calculated in Routine~\ref{rt:tuning} is the minimum value of $\lambda$
with which $\hat{\bmbeta}=\mathbf{0}$ for $\bfy^{(t)}$. Therefore,
under the trial distribution $\pi_r(\bullet;\mathbf{0},(\sigma^2)^{\dag},\lambda^{\dag})$,
$P(\hat{\bmbeta}=\mathbf{0})$ is around 25\% and there is a 75\% of chance
for $\hat{\bmbeta}$ to have some nonzero components. This often results in a good balance between the
dominating region of the target distribution ($\hat{\bmbeta}=\mathbf{0}$) and the region of interest $\Omega_r^*$
for p-value calculation \eqref{eq:pvalue}.
For all numerical examples in this article, we choose $M^{\dag}=5$ and $L_{\text{pilot}}=100$.

\subsection{Multiple tests}

Consider multiple linear models with the same set of predictors,
\begin{equation}\label{eq:multilinear}
\bfy_k=\mathbf{X} \bm{\beta}_k + \bm{\varepsilon}_k, \;\; k=1,\ldots,m,
\end{equation}
where $\bfy_k\in \mathbb{R}^n$, $\bm{\beta}_k \in \mathbb{R}^p$, and $\bm{\varepsilon}_k \sim \mathcal{N}_n(\mathbf{0}, \sigma_k^2 \mathbf{I}_n)$.
After proper rescaling of $\bfy_k$ and $\bm{\beta}_k$, we may assume that all $\sigma_k^2$ are identical,
i.e., $\sigma_k^2=\sigma^2$.
Suppose we are interested in testing against $m$ null hypotheses $\mathcal{H}_k:$
$\bm{\beta}_k=\mathbf{0}$ and $\sigma^2=\sigma_0^2$, given Lasso-type estimates $\hat{\bmbeta}_k^*$ with $\lambda=\lambda_k^*$
for $k=1,\ldots,m$. There are $m$ p-values to calculate,
\begin{equation}\label{eq:multipvalue}
q^*_k=P(|T(\hat{\bmbeta})|\geq T_k^*; \mathcal{H}_k, \lambda_k^*),
\end{equation}
where $T_k^*=|T(\hat{\bmbeta}_k^*)|$
for $k=1,\ldots,m$. This problem occurs in various genomics applications. To give an example, $\bfy_k$
may be the expression level of gene $k$ and $\mathbf{X}$ the expression levels of
$p$ transcription factors across $n$ individuals. The transcription factors may
potentially regulate the expression of a gene
through the linear model \eqref{eq:multilinear}. Rejection of $\mathcal{H}_k$ indicates that gene $k$
is regulated by at least one of the $p$ transcription factors.

To estimate all $q^*_k$, we only need to draw $(\mathbf{b}_A,\mathbf{s}_I,A)^{(t)}$, $t=1,\ldots,L$, from one trial distribution
$\pi_r(\bullet;\mathbf{0},(\sigma^2)^{\dag},\lambda^{\dag})$, in which $(\sigma^2)^{\dag}=M^{\dag} \sigma_0^2$ and
$\lambda^{\dag}$ is obtained by applying the same tuning procedure (Routine~\ref{rt:tuning}) once.
Then we calculate the importance weights by \eqref{eq:ISweight} for all target distributions,
$\{w((\mathbf{b}_A,\mathbf{s}_I,A)^{(t)};\sigma_0^2,\lambda_k^*)\}_{1\leq t \leq L}$, $k=1,\ldots,m$,
and construct estimates for all $q^*_k$ by \eqref{eq:ISest}.

\begin{remark}\label{rm:pathpvalue}
Alternatively, one may apply the Lars algorithm in the direct sampler to draw from the sampling
distribution of $\hat{\bmbeta}$ given $\bm{\beta}=\mathbf{0}$ and $\sigma^2=\sigma_0^2$ for all $\lambda_k^*$,
as the Lars algorithm provides the whole solution path. The computing time of both methods
is dominated by drawing samples and thus is comparable.
However, when $q^*_k$ is small, the IS method will be
orders of magnitude more efficient than direct sampling, and when $q^*_k$ is not too small, the accuracy
of the two methods is on the same order. We will see this in the numerical examples. In addition, we do not have
to use the Lars algorithm to draw from the trial distribution since there is no need to compute the solution path
for importance sampling. One thus has the freedom to choose other algorithms,
such as coordinate descent \citep{Friedman07,WuLange08},
which may be more efficient when both $n$ and $p$ are large.
\end{remark}

\subsection{Numerical examples}

We first simulated two datasets to demonstrate the effectiveness in p-value calculation
by the IS method for individual tests.
Each row of $\mathbf{X}$ was generated from
$\mathcal{N}_p(\mathbf{0},\bm{\Sigma}_{\mathbf{X}})$, where the diagonal and the off-diagonal elements of $\bm{\Sigma}_{\mathbf{X}}$
are 1 and $0.05$, respectively. Given the predictors $\mathbf{X}$, the response vector $\bfy$ was drawn
from $\mathcal{N}_n(\mathbf{X} \bm{\beta}_0, \sigma_0^2 \mathbf{I}_n)$. We set all weights $w_j=1$.
Table~\ref{tab:dataindividual} reports the
values of $n$, $p$, $\sigma_0^2$, and $\bm{\beta}_0$ for the two datasets.
We applied the Lars algorithm on the two datasets and
chose $\lambda^*$ as the first $\lambda$ along the solution path such that
the Lasso estimate $\hat{\bmbeta}^*$ gave the correct number of active coefficients (Table~\ref{tab:dataindividual}).
It turned out that $\hat{\bmbeta}^*$ only included one true active coefficient for both datasets.
Let $A^*=\text{supp}(\hat{\bmbeta}^*)$ be the active set of $\hat{\bmbeta}^*$.
We designed the following test statistics, $T_1=\|\hat{\bmbeta}\|_1$, $T_2=\|\hat{\bmbeta}\|_{\infty}$,
and $\tilde{T}_j=|\hat{\beta}_j|$ for $j \in A^*$, and aimed to calculate p-values under the null hypothesis
$\mathcal{H}_0: \bm{\beta}=\mathbf{0}$ and $\sigma^2=\sigma_0^2$.

\begin{table}[ht]
\centering
\caption{Simulated datasets for individual tests \label{tab:dataindividual}}
\vspace{0.05in}
   \begin{tabular}{c|ccclcc}
   \hline
     Dataset	& $n$ & $p$ & $\sigma_0^2$ & $\bm{\beta}_0$ & $\lambda^*$  & $\lambda^{\dag}$\\
     \hline
     E & 5 & 10 & 1/4 & $(2,-2, 0, \ldots, 0)$ & 1.65 & 0.60\\
     F & 10 & 20 & 1/4 & $(1,1,-1,-1, 0, \ldots,0)$ & 0.315 & 0.57\\
    \hline
  \end{tabular}
\end{table}

We chose $(\sigma^2)^{\dag}=5 \sigma_0^2$ and used Routine~\ref{rt:tuning} to choose $\lambda^{\dag}$
for the trial distributions. The values of $\lambda^{\dag}$ for the two datasets are given in Table~\ref{tab:dataindividual}.
When $(\sigma^2)^{\dag}$ is sufficiently large for a dataset, the $\lambda^{\dag}$ tuned by Routine~\ref{rt:tuning}
can be greater then $\lambda^*$ (dataset F).
The IS method (Routine~\ref{rt:IS}) was applied with $L=5,000$ to estimate p-values for all the above tests.
This estimation procedure was repeated 10 times independently to obtain the standard deviation of an
estimated p-value. We quantify the efficiency of an estimated p-value, $\hat{q}$, by its coefficient of variation
$\text{cv}(\hat{q})=\text{SD}(\hat{q})/\mathbb{E}(\hat{q})$,
where the standard deviation and the mean are calculated across multiple runs. Table~\ref{tab:indpval}
summarizes the results, where $A^*=\{2,3\}$ for dataset E and $A^*=\{4,9,15,17\}$ for dataset F.
One sees that the IS estimates were very accurate: Even for a tail probability as small as $10^{-21}$,
the coefficient of variation was less than or around 2. To benchmark the performance, we approximated
the coefficient of variation of the estimate $\hat{q}^{\text{(DS)}}$ constructed by direct sampling from the target distribution,
$\text{cv}(\hat{q}^{\text{(DS)}})=\sqrt{(1-\bar{q})/(L\bar{q}})$ with $\bar{q}=\mathbb{E}(\hat{q}^{\text{(IS)}})$.
As reported in the table, for estimating an extremely small p-value, $\text{cv}(\hat{q}^{\text{(DS)}})$ can
be orders of magnitude greater than that of an IS estimate, and for a moderate p-value (around $10^{-2}$),
the two methods showed comparable performance.

\begin{table}[t]
\centering
\caption{Estimation of p-values for datasets E and F\label{tab:indpval}}
\vspace{0.05in}
   \begin{tabular}{cc|cccc}
   \hline
  &  & $\mathbb{E}(\hat{q}^{\text{(IS)}})$ & $\text{SD}(\hat{q}^{\text{(IS)}})$ & $\text{cv}(\hat{q}^{\text{(IS)}})$ & $\text{cv}(\hat{q}^{\text{(DS)}})$\\
     \hline
 & $T_1$ & $3.7 \times 10^{-19}$ & $8.7\times 10^{-19}$ & 2.37 & $2.32\times 10^7$ \\
E & $T_2$ & $5.7 \times 10^{-15}$ & $6.2\times 10^{-19}$ & 1.09 & $1.88\times 10^5$ \\
& $\tilde{T}_2$ & $6.3 \times 10^{-21}$ & $1.1\times 10^{-20}$ & 1.81 & $1.79\times 10^8$ \\
& $\tilde{T}_3$ & $8.4 \times 10^{-15}$ & $9.8\times 10^{-15}$ & 1.16 & $1.54\times 10^5$ \\
   \hline
& $T_1$ & $1.5 \times 10^{-6}$ & $4.5\times 10^{-7}$ & 0.30 & 11.5 \\
& $T_2$ & $1.2 \times 10^{-3}$ & $1.2\times 10^{-4}$ & 0.11 & 0.41 \\
F & $\tilde{T}_4$ & $5.7 \times 10^{-5}$ & $2.2\times 10^{-5}$ & 0.38 & 1.87 \\
& $\tilde{T}_9$ & $1.1 \times 10^{-2}$ & $1.3\times 10^{-3}$ & 0.12 & 0.14 \\
& $\tilde{T}_{15}$ & $2.5 \times 10^{-2}$ & $1.9\times 10^{-3}$ & 0.08 & 0.09 \\
& $\tilde{T}_{17}$ & $4.8 \times 10^{-5}$ & $1.5\times 10^{-5}$ & 0.31 & 2.04\\
\hline
  \end{tabular}
\end{table}

Next, we simulated $m=50$ datasets to test our p-value calculation for the multiple testing problem.
We used the design matrix $\mathbf{X}$ in dataset F and $\sigma_0^2=1/4$. The response vector $\bfy_k$
was drawn from $\mathcal{N}_n(\mathbf{X} \bm{\beta}_k, \sigma_0^2 \mathbf{I}_n)$, where the true coefficient vector $\bm{\beta}_k$
is given in Table~\ref{tab:datamultiple} for $k=1,\ldots,50$.
For 10 datasets, $\bm{\beta}_k=\mathbf{0}$ and the null hypothesis is true. For 20 datasets, there are
two large coefficients, which represents the case that the true model is sparse. The other 20 datasets
mimic the scenario in which the true model has many relatively small coefficients.
We chose $\lambda_k^*$ as the first $\lambda$
that gave two active coefficients along the solution path
and used $T(\hat{\bmbeta})=\|\hat{\bmbeta}\|_1$ as the test statistic. Summaries of $\lambda_k^*$ and $T_k^*=\|\hat{\bmbeta}_k^*\|_1$
for the 50 datasets are provided in Table~\ref{tab:datamultiple} as well, from which we see that these datasets
cover a wide range of $\lambda_k^*$ and $T_k^*$.
We chose $(\sigma^2)^{\dag}=5\sigma_0^2$. As seen from Routine~\ref{rt:tuning}, for identical $\mathbf{X}$
and $(\sigma^2)^{\dag}$, the tuning procedure is the same. Therefore, we simply set $\lambda^\dag=0.57$,
the value we used for dataset F (Table~\ref{tab:dataindividual}).

\begin{table}[ht]
\centering
\caption{Simulated datasets for multiple tests \label{tab:datamultiple}}
\vspace{0.05in}
   \begin{tabular}{c|lccc}
   \hline
     Dataset	& $\bm{\beta}_k$ & range of $\lambda_k^*$  & range of $T_k^*$\\
     \hline
     1-10 & $(0, \ldots, 0)$ & (0.16, 0.34) & (0.08, 0.23)\\
     11-30 & $(2,-2, 0, \ldots, 0)$ & (0.88, 1.31) & (0.27, 0.84)\\
     31-50 & $(1/4,\ldots,1/4)$ & (0.70, 1.13) & (0.04, 0.51)\\
    \hline
  \end{tabular}
\end{table}

We simulated $L=5,000$ samples from the trial distribution and estimated the p-values for all the 50 datasets.
This procedure was repeated 10 times independently. The average over 10 runs of the estimated p-value,
$\mathbb{E}(\hat{q}_k^{\text{(IS)}})$, is shown in Figure~\ref{fig:multipval}(a) for $k=1,\ldots,50$.
As expected, most of the p-values for the first 10 datasets
were not significant, while those for the other 40 datasets ranged from $10^{-4}$ to $10^{-30}$,
which confirms that $T(\hat{\bmbeta})=\|\hat{\bmbeta}\|_1$ is a reasonable test statistic.
Again, we see that even for p-values on the order of $10^{-30}$,
the coefficient of variation of an IS estimate was at most around 3 (Figure~\ref{fig:multipval}(b)).
This provides huge gain in accuracy compared to direct sampling. Figure~\ref{fig:multipval}(c) plots
$\log_{10}[\text{cv}(\hat{q}_k^{\text{(DS)}})/\text{cv}(\hat{q}_k^{\text{(IS)}})]$ for the 50 datasets, where
$\text{cv}(\hat{q}_k^{\text{(DS)}})$ was approximated in the same way as in the previous example.
It is comforting to see that while the IS estimates $\hat{q}_k^{\text{(IS)}}$
showed huge improvement over the DS estimates in estimating a tail probability, they were only slightly
worse than the DS estimates for an insignificant p-value. For the first 10 datasets, the coefficient of variation
of $\hat{q}_k^{\text{(IS)}}$ was at most 7.9 times that of $\hat{q}_k^{\text{(DS)}}$. For majority of the other 40 datasets,
the ratio of $\text{cv}(\hat{q}_k^{\text{(DS)}})$ over $\text{cv}(\hat{q}_k^{\text{(IS)}})$ was between
$100$ and $10^{10}$.

\begin{figure}[t]
\centering
   \includegraphics[width=0.55\linewidth,trim=0in 0in 0in 0in,clip]{multipval.eps} \\
   \caption{Estimation of p-values for 50 simulated datasets by IS with a single trial distribution:
    (a) $\log_{10}\mathbb{E}(\hat{q}_k^{\text{(IS)}})$,
   (b) $\text{cv}(\hat{q}_k^{\text{(IS)}})$, and (c) $\log_{10}[\text{cv}(\hat{q}_k^{\text{(DS)}})/\text{cv}(\hat{q}_k^{\text{(IS)}})]$ for $k=1,\ldots,50$. Each bar in a plot gives the result for one dataset.
    \label{fig:multipval}}
\end{figure}

\section{Estimating sampling distributions}\label{sec:errorbounds}

For a vector $\mathbf{v}=(v_j)$, denote its active set by $\mathcal{A}(\mathbf{v})=\{j:v_j\ne 0\}$. Let $\bm{\beta}_0=(\beta_{0j})_{1:p}$ be the true coefficient vector and $A_0=\mathcal{A}(\bm{\beta}_0)$. Consider model \eqref{eq:linearmodel} with $\bm{\beta}=\bm{\beta}_0$ and $\bm{\varepsilon} \sim \mathcal{N}_n(\mathbf{0},\sigma^2\mathbf{I}_n)$. The penalized loss in \eqref{eq:lassoloss} is, up to an additive constant,
\begin{align}
 &\, \frac{1}{2} \|\bfy-\mathbf{X} \bm{\beta} \|_2^2 -\frac{1}{2} \|\bfy-\mathbf{X} \bm{\beta}_0 \|_2^2
 + n \lambda \sum_{j=1}^p w_j (| \beta_{j}|-|\beta_{0j}|)  \nonumber\\
= &\, \frac{n}{2r_n^2} \bm{\delta}^{\mathsf{T}} \mathbf{C} \bm{\delta} - \frac{n}{r_n}\bm{\delta}^{\mathsf{T}} \mathbf{U}
+ n\lambda \sum_{j \in A_0} w_j \left(| \beta_{0j}+ r_n^{-1}\delta_{j}|- | \beta_{0j}|\right)
+ \frac{n\lambda}{r_n} \sum_{j \notin A_0} w_j |\delta_{j} | \nonumber\\
\defi &\, V(\bm{\delta}; \bm{\beta}_0, \mathbf{U}), \label{eq:defVn}
\end{align}
where $r_n>0$, $\bm{\delta} = (\delta_j)=r_n (\bm{\beta} - \bm{\beta}_0)$, and $\mathbf{U}=\mathbf{X}^{\mathsf{T}}\bm{\varepsilon}/n$.
Assuming $V$ has a unique minimizer,
\begin{equation}\label{eq:defdelta}
\arg\min_{\bm{\delta}} V(\bm{\delta};\bm{\beta}_0,\mathbf{U})=r_n(\hat{\bmbeta}-\bm{\beta}_0)\defi \hat{\bm{\delta}},
\end{equation}
where $\hat{\bmbeta}$ is the Lasso-type estimator that minimizes \eqref{eq:lassoloss}. Suppose that $\check{\bmbeta}$ and $\hat{\sigma}$ are estimators of $\bm{\beta}_0$ and $\sigma$, following their respective sampling distributions. Let
\begin{equation}\label{eq:defdeltastar}
\bmdelta^*=r_n({\bmbeta}^*-\check{\bmbeta})=\arg\min_{\bm{\delta}} V(\bm{\delta};\check{\bmbeta},\mathbf{U}^*),
\end{equation}
where $\mathbf{U}^*=\mathbf{X}^{\mathsf{T}}\bm{\varepsilon}^*/n$ and $\bm{\varepsilon}^*\mid \hat{\sigma} \sim \mathcal{N}_n(\mathbf{0},\hat{\sigma}^2\mathbf{I}_n)$. For random vectors $\mathbf{Z}_1$ and $\mathbf{Z}_2$, we use $\nu[\mathbf{Z}_1]$ to denote the distribution of $\mathbf{Z}_1$ and $\nu[\mathbf{Z}_1\mid\mathbf{Z}_2]$ to denote the conditional distribution of $\mathbf{Z}_1$ given $\mathbf{Z}_2$. In general, $\nu[\mathbf{Z}_1\mid\mathbf{Z}_2]$ is a random probability measure. The goal of this section is to derive nonasymptotic bounds on the difference between $\nu[\hat{\bm{\delta}}]$ and $\nu[\bmdelta^*\mid \check{\bmbeta},\hat{\sigma}]$, with all proofs relegated to Section~\ref{sec:proofs}. Our results provide theoretical justifications for any method that estimates the uncertainty in $\hat{\bmbeta}$ via simulation from $\nu[\bmdelta^*\mid \check{\bmbeta},\hat{\sigma}]$. In particular, for estimator augmentation, we draw $({\bmbeta}^*,\mathbf{S}^*)|\check{\bmbeta}, \hat{\sigma}^2$ from an estimated sampling distribution whose density is given by \eqref{eq:jointdensitynormal} or \eqref{eq:densitynormalhigh} with $\bm{\beta}=\check{\bmbeta}$ and $\sigma^2=\hat{\sigma}^2$.

\subsection{Relevant existing results}

We compile relevant published results here. The following restricted eigenvalue (RE) assumption is a special case of the one used in \cite{Lounici11}, which extends the original definition in \cite{Bickel09}.

\begin{assumption}[RE$(m,c_0)$]
For some positive integer $m \leq p$ and a positive number $c_0$, the following condition holds:
\begin{equation*}
\kappa(m,c_0) \defi \min_{|A| \leq m} \min_{\bm{\delta} \ne \mathbf{0}}
\left\{\frac{\|\mathbf{X} \bm{\delta}\|_2}{\sqrt{n} \|\bm{\delta}_{A} \|_2}:
\sum_{j\in A^c} w_j |\delta_{j}|  \leq c_0 \sum_{j\in A} w_j | \delta_{j}|\right\} >0.
\end{equation*}
\end{assumption}

Lemma~\ref{lm:lassobounds} is from Theorem 3.1 and its proof in \cite{Lounici11}, regarding the Lasso as a special case of the group Lasso \citep{YuanLin06} with the size of every group being one.
\begin{lemma}\label{lm:lassobounds}
Consider the model \eqref{eq:linearmodel} with $\bm{\beta}=\bm{\beta}_0$ and $\bm{\varepsilon} \sim \mathcal{N}_n(\mathbf{0},\sigma^2\mathbf{I}_n)$, $\sigma^2>0$, and let $p\geq 2$, $n \geq 1$. Suppose that all the diagonal elements of $\mathbf{C}$ are $1$ and $\mathbf{U}=(U_j)_{1:p}=\frac{1}{n}\mathbf{X}^{\mathsf{T}}\bm{\varepsilon}$. Let $q_0= |A_0|\leq q$, where $1\leq q \leq p$, and Assumption RE$(q,3)$ be satisfied. Choose $u>1$ and
\begin{align}\label{eq:condlambda}
\lambda \geq  \lambda_0 \defi\frac{2\sigma}{w_{\min}} \sqrt{(2+5u\log p)/n},
\end{align}
where $w_{\min}=\inf\{w_j: j=1,\ldots,p\}$.
Then on the event $\mathcal{E}=\cap_{j=1}^p\{|U_j|\leq w_j \lambda/2\}$, which happens with probability at least $1-2p^{1-u}$, for any minimizer $\hat{\bmbeta}$ of \eqref{eq:lassoloss} we have
\begin{align}
|\mathcal{A}(\hat{\bmbeta})| & \leq \frac{64\phi_{\max}}{\kappa^2(q,3)} \sum_{A_0} \frac{w_j^2}{w_{\min}^2},  \label{eq:boundsize}\\
\|\hat{\bmbeta}-\bm{\beta}_0 \|_1 & \leq \frac{16\lambda}{\kappa^2(q,3)} \sum_{A_0} \frac{w_j^2}{w_{\min}},
\label{eq:boundL1}
\end{align}
where $\phi_{\max}$ is the maximum eigenvalue of $\mathbf{C}$. If Assumption RE$(2q,3)$ is satisfied, then on the same event,
\begin{align}
\|\hat{\bmbeta}-\bm{\beta}_0 \|_2 & \leq \frac{4\sqrt{10}}{\kappa^2(2q,3)} \frac{\lambda \sum_{A_0}w_j^2}{w_{\min}\sqrt{q}}\defi \tau.
\label{eq:boundL2}
\end{align}
\end{lemma}

As an immediate consequence of this lemma, the Lasso-type estimator has the screening property assuming a suitable beta-min condition.
\begin{lemma}\label{lm:screen}
Let the assumptions in Lemma~\ref{lm:lassobounds} be satisfied. If $\inf_{A_0} | \beta_{0j}|>\tau$,
then on the event $\mathcal{E}$, $A_0\subseteq \mathcal{A}(\hat{\bmbeta})$ for any $\hat{\bmbeta}$. If $\inf_{A_0} | \beta_{0j}| > 2\tau$,
then on the same event,
\begin{equation*}
\left\{j: |\hat{\beta}_j|>\tau\right\} = A_0.
\end{equation*}
\end{lemma}

\subsection{Known variance}\label{sec:errorboundsknownvar}

In this subsection we assume that the noise variance $\sigma^2$ is known and fix $\hat{\sigma}=\sigma$ in the definition of $\bmdelta^*$ \eqref{eq:defdeltastar}. We first regard $\check{\bmbeta}=(\check{\beta}_j)_{1:p}$ as a fixed vector and find conditions which are sufficient for the distribution of $\bmdelta^*$ to be close to that of $\hat{\bm{\delta}}$. Then we construct an estimator that satisfies these conditions with high probability. Let
\begin{equation}\label{eq:difftdbeta}
\eta  \defi \sup_{j \in A_0} \frac{|\check{\beta}_{j}-\beta_{0j}|}{| \beta_{0j}|}.
\end{equation}

\begin{lemma}\label{lm:fixknownvar}
Assume that the columns of $\mathbf{X}$ are in general position. Fix $\hat{\sigma}=\sigma$. Suppose that $\check{\bmbeta}$ is a fixed vector in $\mathbb{R}^p$ so that (i) $\check{\beta}_{j}=0$ for all $j \notin A_0$ and (ii) $\eta\in[0,1)$ as defined in \eqref{eq:difftdbeta}. Let $M_1>0$ and assume
\begin{align}\label{eq:assbeta0j}
\inf_{j\in A_0}|\beta_{0j}|>  \frac{M_1}{r_n(1-\eta)}.
\end{align}
Then we have
\begin{equation*}
\nu[\bmdelta^*\mid \|\bmdelta^* \|_{\infty}< M_1]=\nu[\hat{\bm{\delta}}\mid \|\hat{\bm{\delta}}\|_{\infty}<M_1].
\end{equation*}
\end{lemma}

One possible way to construct $\check{\bmbeta}$ that satisfies conditions (i) and (ii) in Lemma~\ref{lm:fixknownvar} is to threshold the Lasso-type estimator $\hat{\bmbeta}$ by a constant $b_{\text{th}}>0$, i.e.,
\begin{equation}\label{eq:thresholddef}
\check{\beta}_{j}=\hat{\beta}_{j}\bm{1}(|\hat{\beta}_{j}|>b_{\text{th}}), \quad j=1,\ldots,p.
\end{equation}

\begin{theorem}\label{thm:threshold}
Let the assumptions in Lemma~\ref{lm:lassobounds} be satisfied and assume that the columns of $\mathbf{X}$ are in general position. Choose $r_n$ such that $\|\hat{\bm{\delta}}\|_{\infty} < M_1$ with probability at least $1-\alpha_1$. Fix $\hat{\sigma}=\sigma$, define $\check{\bmbeta}$ by \eqref{eq:thresholddef} with $b_{\text{th}}=\tau$ \eqref{eq:boundL2}, and assume
\begin{align}\label{eq:betamin}
\inf_{j\in A_0}|\beta_{0j}|>  \max\left\{2\tau, \frac{M_1}{r_n(1-\eta)}\right\}.
\end{align}
Then with probability at least $1-2p^{1-u}$, we have
\begin{equation}
\sup_{B\in \mathscr{R}^p} |P(\bmdelta^* \in B \mid \check{\bmbeta}) - P(\hat{\bm{\delta}} \in B)| \leq 2\alpha_1, \label{eq:supPrdiff}
\end{equation}
where $\mathscr{R}^p$ is the $\sigma$-field of $p$-dimensional Borel sets.
\end{theorem}

\begin{remark}
Depending on the estimator $\check{\bmbeta}$, the conditional probability $P(\bmdelta^* \in B \mid \check{\bmbeta})$ is a random variable. The probability $1-2p^{1-u}$ is with respect to the sampling distribution of $\bfy$. This theorem gives an explicit nonasymptotic bound \eqref{eq:supPrdiff} on the difference between the two probability measures, $\nu[\bmdelta^* \mid \check{\bmbeta}]$ and $\nu[\hat{\bm{\delta}}]$, and justifies simulation from the estimated sampling distribution with $\bm{\beta}=\check{\bmbeta}$.
\end{remark}

\begin{remark}\label{rmk:asymknownvar}
Consider the asymptotic implications of Theorem~\ref{thm:threshold} by allowing $q_0,p\to \infty$ as $n\to \infty$. Suppose that $w_{\min}$ and $w_{\max}=\sup_j w_j$ stay bounded away from 0 and $\infty$. For the Lasso, $w_{\min}=w_{\max}=1$. Choose $r_n$ so that $\|\hat{\bm{\delta}}\|_{\infty}\asymp_P 1$, i.e., the convergence rate of $\|\hat{\bmbeta}-\bm{\beta}_0\|_{\infty}$ is $1/r_n$. For any $\alpha_1>0$, there is $M_1<\infty$ so that $P(\|\hat{\bm{\delta}}\|_{\infty}<M_1)\geq 1-\alpha_1$. Let $q=q_0$ in Lemma~\ref{lm:lassobounds} and assume that $\liminf_n\kappa(m,3)>0$ when $m=O(q_0)=o(n)$. With a suitable choice of $\lambda \asymp \sqrt{(\log p)/n}$,
\begin{equation*}
\tau \asymp \sqrt{q_0} \lambda \asymp \sqrt{q_0 (\log p)/n}.
\end{equation*}
Because $1 \asymp_P r_n\|\hat{\bmbeta}-\bm{\beta}_0\|_{\infty}=O_P(r_n\tau)$, we have $r_n^{-1}=O(\tau)$. Thus, the order of $r_n$ satisfies
\begin{equation}\label{eq:orderofr_n}
r_n^{-1}=O(\sqrt{q_0 (\log p)/n}) \text{ and } r_n = O(\sqrt{n}).
\end{equation}
Consequently, a sufficient condition for \eqref{eq:betamin} is
\begin{equation}\label{eq:orderofbetamin}
\inf_{A_0}|\beta_{0j}|\gg  \sqrt{q_0 (\log p)/n}
\end{equation}
with the order of $b_{\text{th}}$ in between. Theorem~\ref{thm:threshold} then implies that \eqref{eq:supPrdiff} holds with probability at least $1-2p^{1-u}\to 1$ as $p\to \infty$. Choosing $\alpha_1$ arbitrarily close to zero, this demonstrates that
\begin{equation}\label{eq:supPrdifftinP}
\sup_{B\in \mathscr{R}^p} |P(\bmdelta^* \in B \mid \check{\bmbeta}) - P(\hat{\bm{\delta}}\in B )|\toinP 0.
\end{equation}
As $\inf_{A_0}|\beta_{0j}|$ may decay to zero at a rate slower than $\sqrt{{q_0 (\log p)}/{n}}\to 0$, Theorem~\ref{thm:threshold} applies in the high-dimensional setting ($p\gg n\to\infty$).
\end{remark}

\subsection{Unknown variance}\label{sec:theoryunknownvar}

When $\sigma^2$ is unknown, recall that $\bm{\varepsilon}^*\mid \hat{\sigma} \sim \mathcal{N}_n(\mathbf{0},\hat{\sigma}^2\mathbf{I}_n)$ and $\mathbf{U}^*=\mathbf{X}^{\mathsf{T}}\bm{\varepsilon}^*/n$. Denote the components of $\mathbf{U}^*$ by $U_j^*$, $j=1,\ldots,p$. Define
\begin{align}\label{eq:defdeltastar0}
\bmdelta^*_0=\arg\min_{\bm{\delta}} V(\bm{\delta};\check{\bmbeta},(\sigma/\hat{\sigma})\mathbf{U}^*),
\end{align}
whose distribution does not depend on $\hat{\sigma}$ and is identical to that of $\bmdelta^*$ when $\hat{\sigma}$ is fixed to $\sigma$. Thus, results in Section~\ref{sec:errorboundsknownvar} show that $\nu[\bmdelta^*_0\mid \check{\bmbeta}]$ is close to $\nu[\hat{\bm{\delta}}]$. We will further bound the difference, $\bmdelta^*-\bmdelta^*_0$, in this subsection. In other words, $\bmdelta^*_0$ serves as an intermediate variable between $\hat{\bm{\delta}}$ and $\bmdelta^*$ to help quantify the difference between their distributions.

\begin{lemma}\label{lm:fixunknownvar}
Assume that the columns of $\mathbf{X}$ are in general position. Let $\check{\bmbeta}$ and $\hat{\sigma}$ be fixed such that $|\mathcal{A}(\check{\bmbeta})|\leq q$ with $1\leq q \leq p$ and
\begin{equation}\label{eq:diffhsigma}
\max\{\left|{\hat{\sigma}}/{\sigma}-1\right|, \left|{\sigma}/{\hat{\sigma}}-1\right|\} \leq \zeta \in [0,1).
\end{equation}
If Assumption RE$(q,3)$ is satisfied, then on the event $\cap_{j=1}^p\{|U_j^*| \leq (1-\zeta)w_j\lambda/2\}$,
\begin{equation}\label{eq:boundquadratic}
\frac{1}{n}\|\mathbf{X} (\bmdelta^*-\bmdelta^*_0) \|^2_2 \leq \frac{32 w_{\max}r_n^2\lambda^2 \zeta}{w_{\min}\kappa^2(q,3)}  \sum_{A_0} w_j^2.
\end{equation}
\end{lemma}

Based on \eqref{eq:boundquadratic}, we can obtain an upper bound on $\|\bmdelta^*-\bmdelta^*_0\|_2$ via the restricted eigenvalues of $\mathbf{C}$. Define the minimum restricted eigenvalue of $\mathbf{C}$ for an integer $m\leq p$ by
\begin{equation}\label{eq:remin}
\phi_{\min}(m)=\min_{1\leq |\mathcal{A}(\mathbf{v})|\leq m} \frac{\mathbf{v}^{\mathsf{T}}\mathbf{C} \mathbf{v}}{\|\mathbf{v}\|_2^2},
\end{equation}
and let
\begin{equation*}
M_2=\frac{128\phi_{\max}}{\kappa^2(q,3)} \sum_{A_0} \frac{w_j^2}{w_{\min}^2},
\end{equation*}
which is twice the upper bound on $|\mathcal{A}(\hat{\bmbeta})|$ in \eqref{eq:boundsize}.

\begin{theorem}\label{thm:thresholdsigmaest}
Let the assumptions in Theorem~\ref{thm:threshold} be satisfied but without fixing $\hat{\sigma}$ to $\sigma$. In addition, assume that \eqref{eq:diffhsigma} holds with probability at least $1-\alpha_2$ and $\phi_{\min}(M_2)>0$. Choose
\begin{equation}\label{eq:condlambdasigmaest}
\lambda \geq \frac{1+\zeta}{1-\zeta} \lambda_0,
\end{equation}
where $\lambda_0$ is defined in \eqref{eq:condlambda}.
Then with probability at least $1-(\alpha_2+2p^{1-u})$, we have
\begin{equation}\label{eq:supPrdiffdelta0}
\sup_{B\in \mathscr{R}^p} |P(\bmdelta^*_0 \in B \mid \check{\bmbeta}) - P(\hat{\bm{\delta}} \in B)| \leq 2\alpha_1
\end{equation}
and
\begin{equation}\label{eq:L2Deltabound}
P\left\{\left.\|\bmdelta^*-\bmdelta^*_0 \|^2_2 \leq \frac{32 w_{\max}r_n^2\lambda^2 \zeta \sum_{A_0} w_j^2}{w_{\min}\kappa^2(q,3)\phi_{\min}(M_2)}\right| \check{\bmbeta}, \hat{\sigma}\right\}\geq 1- 2p^{1-u}.
\end{equation}
\end{theorem}

\begin{remark}
Assume that $\max\{\left|{\hat{\sigma}}/{\sigma}-1\right|, \left|{\sigma}/{\hat{\sigma}}-1\right|\} =O_P(\zeta_n)$ with $\zeta_n\to0$. Following the asymptotic framework in Remark~\ref{rmk:asymknownvar}, we can establish by the same reasoning that
\begin{equation}\label{eq:uniformdelta0}
\sup_{B\in \mathscr{R}^p} |P(\bmdelta^*_0 \in B \mid \check{\bmbeta}) - P(\hat{\bm{\delta}}\in B )|\toinP 0.
\end{equation}
Suppose that $\liminf_n\phi_{\min}(m)>0$ for $m=o(n)$ and that $\phi_{\max}$ is bounded from above by a constant (or at least does not diverge too fast). Then $M_2=O(\phi_{\max}q_0)=o(n)$ and $\phi_{\min}(M_2)$ is bounded from below by a positive constant. The upper bound in \eqref{eq:L2Deltabound} becomes $O(r_n^2 \lambda^2 \zeta_n q_0)=O(\zeta_n q_0 \log p)$ as $r_n=O(\sqrt{n})$ \eqref{eq:orderofr_n} and $\lambda^2\asymp \log(p)/n$. If $\zeta_n =o(1/(q_0 \log p))$, then \eqref{eq:L2Deltabound} implies that
\begin{equation*}
P(\|\bmdelta^*-\bmdelta^*_0\|_2>\epsilon\mid\check{\bmbeta},\hat{\sigma})\to 0
\end{equation*}
for any $\epsilon>0$. Combing with \eqref{eq:uniformdelta0} this shows that, with probability tending to one, $\nu[\bmdelta^*_J\mid \check{\bmbeta},\hat{\sigma}]$ converges weakly to $\nu[\hat{\bm{\delta}}_J]$ for any fixed index set $J\subseteq\{1,\ldots,p\}$. If $\sigma$ is estimated by the scaled Lasso \citep{SunZhang12}, one may reach $\zeta_n=O(q_0(\log p)/n)$ under certain conditions by their Theorem 2 and it is then sufficient to have $q_0\log p\ll \sqrt{n}$. If $p$ is fixed, we only need $\zeta_n =o(1)$ and hence any consistent estimator of $\sigma$ will be sufficient. In this case, $r_n\asymp \sqrt{n}$ \eqref{eq:orderofr_n} and with probability tending to one, $\nu[\bmdelta^*\mid \check{\bmbeta},\hat{\sigma}]$ converges weakly to $\nu[\hat{\bm{\delta}}]$.
\end{remark}

The key assumptions on the underlying model are the RE assumption on the Gram matrix $\mathbf{C}$ and beta-min and sparsity assumptions on the true coefficients $\bm{\beta}_0$, which are comparable to those in Lemma~\ref{lm:lassobounds} and Lemma~\ref{lm:screen}. There is an extra assumption that $\phi_{\min}(M_2)>0$ in Theorem~\ref{thm:thresholdsigmaest}, which is again imposed on the restricted eigenvalues of $\mathbf{C}$. If $\mathbf{X}$ is drawn from a continuous distribution on $\mathbb{R}^{n\times p}$, then with probability one $\phi_{\min}(m)>0$ for any $m\leq n$. It should be noted that we do not assume the irrepresentable condition \citep{ZhaoYu06,Meinshausen06,Zou06}, which is much stronger than the assumptions on the restricted eigenvalues of $\mathbf{C}$.

We compare our results to residual bootstrap for approximating the sampling distribution of the Lasso. To be precise, a residual bootstrap is equivalent to Routine~\ref{rt:directsample} with $\bm{\beta}$ estimated by $\check{\bmbeta}$ and $\bm{\varepsilon}^{(t)}$ drawn by resampling residuals. Assuming $p$ is fixed, \cite{KnightFu00} argue that the residual bootstrap may be consistent if $\check{\bmbeta}$ is model selection consistent and \cite{Chatterjee11} establish such fixed-dimensional consistency when $\check{\bmbeta}$ is constructed by thresholding the Lasso, in the same spirit as \eqref{eq:thresholddef}. As discussed above, Theorem~\ref{thm:thresholdsigmaest} applied to a fixed $p$ is clearly in line with these previous works. However, our results are much more general by providing explicit nonasymptotic bounds that imply consistency when $p\gg n\to \infty$. More recently, \cite{Chatterjee13} have shown that the residual bootstrap is consistent for the adaptive Lasso \citep{Zou06} when $p>n\to \infty$ under a number of conditions. A fundamental difference is that the weights $w_j$ are specified by an initial $\sqrt{n}$-consistent estimator in their work and will not stay bounded as $n\to\infty$. Therefore, their results do not apply to the Lasso. On the contrary, the results in this section are derived assuming the weights $w_j$ are constants without being specified by any initial estimator. In addition, \cite{Chatterjee13} impose in their Theorem 5.1 that $\inf_{A_0}|\beta_{0j}|\geq K$ for some $K\in(0,\infty)$ when $p>n$, which disallows the decay of the magnitudes of nonzero coefficients as $n$ grows. This is considerably stronger than our assumption \eqref{eq:orderofbetamin}.

\section{Generalizations and discussions}\label{sec:generalization}

\subsection{Random design}

We generalize the Monte Carlo methods to a random design,
assuming that $\mathbf{X}$ is drawn from a distribution $f_{\mathbf{X}}$.
The distribution of the augmented estimator $(\hat{\bmbeta}_{\mathcal{A}},\mathbf{S}_{\mathcal{I}},\mathcal{A})$, \eqref{eq:jointdensity}
and \eqref{eq:jointdensityhigh}, becomes a conditional distribution given $\mathbf{X}=\bfx$, written as
$\pi(\mathbf{b}_A,\mathbf{s}_I,A\mid \bfx)$ and ${\pi}_r(\mathbf{b}_A,\mathbf{s}_I,A\mid \bfx)$, respectively.

In the low-dimensional setting,
we may generalize the MLS (Routine~\ref{rt:MHLS}) to draw samples from
${\pi}(\mathbf{b}_A,\mathbf{s}_I,A,\bfx)={\pi}(\mathbf{b}_A,\mathbf{s}_I,A\mid \bfx){f}_{\mathbf{X}}(\bfx)$
and approximate the sampling distribution of $\hat{\bmbeta}$.
It may be difficult to assume or estimate a reliable density for $\mathbf{X}$,
but it is sufficient for the development of an MH sampler under a random design (rdMLS) if we can draw from ${f}_{\mathbf{X}}(\bfx)$.
As seen below, we do not need an
explicit form of $f_{\mathbf{X}}(\bfx)$ for computing the MH ratio \eqref{eq:rMHration}
and thus may draw $\bfx^{\dag}$ by the bootstrap.

\begin{routine}[rdMLS]\label{rt:rMHLS}
Suppose the current sample is $(\mathbf{b}_A,\mathbf{s}_I,A,\bfx)^{(t)}$.
\begin{itemize}
\item[(1)]{Draw $\bfx^{\dag}$ from ${f}_{\mathbf{X}}$, and accept it as $\bfx^{(t+1)}$ with probability
\begin{equation}\label{eq:rMHration}
\min\left\{1, \frac{\pi((\mathbf{b}_A,\mathbf{s}_I,A)^{(t)}\mid \bfx^{\dag})}
{\pi((\mathbf{b}_A,\mathbf{s}_I,A)^{(t)}\mid \bfx^{(t)})} \right\};
\end{equation}
otherwise, set $\bfx^{(t+1)}=\bfx^{(t)}$.}
\item[(2)]{Regarding ${\pi}(\mathbf{b}_A,\mathbf{s}_I,A\mid \bfx^{(t+1)})$ as the target density,
apply one iteration of the MLS (Routine~\ref{rt:MHLS}) to obtain $(\mathbf{b}_A,\mathbf{s}_I,A)^{(t+1)}$.}
\end{itemize}
\end{routine}

Generalization of the IS algorithm (Routine~\ref{rt:IS}) is also straightforward.
Draw $\bfx^{(t)}$ from $f_{\mathbf{X}}$ and draw $(\mathbf{b}_A,\mathbf{s}_I,A)^{(t)}$ from the trial distribution
given $\mathbf{X}=\bfx^{(t)}$. Calculate importance weights by \eqref{eq:ISweight} with $\mathbf{X}=\bfx^{(t)}$,
and apply the same estimation \eqref{eq:ISest}. Again, an explicit expression for $f_{\mathbf{X}}$ is unnecessary.
But bootstrap sampling from $\mathbf{X}$ is not a choice for the high-dimensional setting,
because a bootstrap sample from $\mathbf{X}$ violates Assumption~\ref{as:X}.

\subsection{Model selection consistency}

The distribution of the augmented estimator may help
establish asymptotic properties of a Lasso-type
estimator. Here, we demonstrate this point by studying the model selection consistency of
the Lasso. Our goal is not to establish new asymptotic results, but to provide an intuitive
and geometric understanding of the technical conditions in existing work.
Recall that $\mathcal{A}$ and $A_0$ are the respective active sets of $\hat{\bmbeta}$ and $\bm{\beta}_0$.
Let $q_0=|A_0|$ and $\mathbf{s}_{0}=\mbox{sgn}(\bm{\beta}_{0A_0})$.
Without loss of generality, assume $A_0=\{1,\ldots,q_0\}$ and $I_0=\{q_0+1,\ldots,p\}$.
We allow both $p$ and $q_0$ to grow with $n$.

\begin{definition}[sign consistency \citep{MeinshausenYu09}]\label{def:signconsistency}
We say that $\hat{\bmbeta}$ is sign consistent for $\bm{\beta}_0$ if
\begin{equation}\label{eq:signconsistency}
P(\mathcal{A}=A_0, \mbox{sgn}(\hat{\bmbeta}_{A_0})=\mathbf{s}_0) \to 1, \text{ as } n\to \infty.
\end{equation}
\end{definition}

If $\hat{\bmbeta}$ is unique, the size of its active set $|\mathcal{A}|\leq n$ (Lemma~\ref{lm:rank}),
and thus $\mathbf{D}(\mathcal{A})$ is invertible from \eqref{eq:defD}.
Therefore, the definitions of $\bm{\mu}(A,\mathbf{s}_A;\bm{\beta})$ and $\bm{\Sigma}(A;\sigma^2)$ in \eqref{eq:munorm}
and \eqref{eq:varnorm} are also valid for any $(\mathbf{b}_A,\mathbf{s}_I,A)\in \Omega_r$ when $p>n$.
Rewrite the KKT condition \eqref{eq:lassojoint} as
\begin{equation}\label{eq:defTheta}
\bm{\Theta}=[\mathbf{D}(\mathcal{A})]^{-1} \mathbf{U} + \bm{\mu}(\mathcal{A},\mathbf{S}_{\mathcal{A}};\bm{\beta}_0),
\end{equation}
where $\bm{\Theta}=(\hat{\bmbeta}_{\mathcal{A}},\mathbf{S}_{\mathcal{I}}) \in \mathbb{R}^p$.
Fixing $\mathcal{A}=A_0$ and $\mathbf{S}_{A_0}=\mathbf{s}_0$ in \eqref{eq:defTheta}, we define a random vector
\begin{equation}\label{eq:defZ}
\mathbf{Z}=[\mathbf{D}(A_0)]^{-1} \mathbf{U} + \bm{\mu}(A_0,\mathbf{s}_0;\bm{\beta}_0)
\end{equation}
via an affine map of $\mathbf{U}$.
Note that we always have $\mathbb{E}(\mathbf{U})=\mathbf{0}$ and $\text{Var}(\mathbf{U})=\frac{\sigma^2}{n}\mathbf{C} \geq 0$, regardless of the sizes of $n$ and $p$.
When $p>n$, $\text{Var}(\mathbf{U})$ is semipositive definite,
meaning that components of $\mathbf{U}$ are linearly dependent of each other, since $\mathbf{U}$ only lies in $\text{row}(\mathbf{X})$,
a proper subspace of $\mathbb{R}^p$. Consequently, $\mathbb{E}(\mathbf{Z})=\bm{\mu}(A_0,\mathbf{s}_{0};\bm{\beta}_0)\defi \bm{\mu}^0$ and
$\text{Var}(\mathbf{Z})=\bm{\Sigma}(A_0;\sigma^2) \defi \bm{\Sigma}^0\geq 0$.
Simple calculation from \eqref{eq:munorm} and \eqref{eq:varnorm} gives
\begin{eqnarray}
&\left(\begin{array}{c}
\bm{\mu}^0_{A_0} \\ \bm{\mu}^0_{I_0} \end{array}\right)  =
\left(
\begin{array}{c}
\bm{\beta}_{0A_0}-\lambda\mathbf{C}_{A_0A_0}^{-1}\mathbf{W}_{A_0A_0} \mathbf{s}_{0}  \\
\mathbf{W}_{I_0I_0}^{-1}\mathbf{C}_{I_0A_0}\mathbf{C}_{A_0A_0}^{-1}\mathbf{W}_{A_0A_0}\mathbf{s}_{0}
\end{array}
\right), \label{eq:muA0} \\
&\bm{\Sigma}^0 =
\frac{\sigma^2}{n}\left( \begin{array}{cc}
\mathbf{C}_{A_0A_0}^{-1} &  \mathbf{0} \\
\mathbf{0} & \lambda^{-2} \mathbf{W}_{I_0I_0}^{-1}\mathbf{C}_{I_0\mid A_0} \mathbf{W}_{I_0I_0}^{-1}
 \end{array}\right), \label{eq:sigmaA0}
\end{eqnarray}
where $\mathbf{C}_{I_0\mid A_0} = \mathbf{C}_{I_0I_0}-\mathbf{C}_{I_0A_0}\mathbf{C}_{A_0A_0}^{-1}\mathbf{C}_{A_0I_0}$.

\begin{lemma}\label{lm:signtoZ}
If the columns of $\mathbf{X}$ are in general position and $\bm{\varepsilon}$ is i.i.d. with mean zero and variance $\sigma^2$, then for any $p\geq 1$ and $n \geq 1$,
\begin{equation}\label{eq:normsign}
P(\mathcal{A}=A_0, \mbox{sgn}(\hat{\bmbeta}_{A_0})=\mathbf{s}_0) = P(\mathbf{Z} \in \Omega_{A_0,\mathbf{s}_0}),
\end{equation}
where $\Omega_{A_0,\mathbf{s}_0}$ is defined by \eqref{eq:subsetAs}.
\end{lemma}
\begin{proof}
If $\mathcal{A}=A_0$ and $\mathbf{S}_{A_0}=\mathbf{s}_0$, then $\bm{\Theta}\in\Omega_{A_0,\mathbf{s}_0}$ by definition. In this case, \eqref{eq:defTheta} reduces to \eqref{eq:defZ} and we have $\mathbf{Z}=\bm{\Theta} \in \Omega_{A_0,\mathbf{s}_0}$. Reversely, if $\mathbf{Z} \in \Omega_{A_0,\mathbf{s}_0}$, then $(\mathbf{Z},A_0)$ is a solution to the KKT condition \eqref{eq:defTheta}. By uniqueness, $(\bm{\Theta},\mathcal{A})=(\mathbf{Z},A_0)$ and therefore $\mathcal{A}=A_0$ and $\mathbf{S}_{A_0}=\mathbf{s}_0$.
\end{proof}

Consequently, to establish sign consistency, we only need a set of sufficient conditions for
$P(\mathbf{Z} \in \Omega_{A_0,\mathbf{s}_0})\to 1$:
(C1) $\mbox{sgn}(\bm{\mu}^0_{A_0})= \mathbf{s}_0$.
(C2) $\|\bm{\mu}^0_{I_0} \|_{\infty} \leq c$ for some $c \in (0,1)$.
(C3) Let $Z_j$ and $\mu^0_j$ be the $j^\text{th}$ components of $\mathbf{Z}$ and $\bm{\mu}^0$, respectively.
As $n\to\infty$,
\begin{equation}\label{eq:C3}
P\left(|Z_j - \mu^0_j| < \delta_j, \,\forall j \right)\to 1,
\end{equation}
where $\delta_j=|\mu^0_j|$ for $j\in A_0$ and $\delta_j=1-c$ for $j\in I_0$.

The first two conditions ensure that
$\bm{\mu}^0=\mathbb{E}(\mathbf{Z})$ lies in the interior of $\Omega_{A_0,\mathbf{s}_0}$. The third condition guarantees that $\mathbf{Z}$
always stays in a box centered at $\bm{\mu}^0$, and the box is contained in $\Omega_{A_0,\mathbf{s}_0}$ if (C1) and (C2) hold.
These conditions have a simple and intuitive geometric interpretation illustrated in Figure~\ref{fig:consistent}.
\begin{lemma}\label{lm:signconsistency}
Assume the columns of $\mathbf{X}$ are in general position and $\bm{\varepsilon}$ is i.i.d. with mean zero and variance $\sigma^2$. If conditions
(C1), (C2), and (C3) hold as $n\to\infty$, then $\hat{\bmbeta}$ is sign consistent for $\bm{\beta}_0$, regardless of the relative
size between $p$ and $n$.
\end{lemma}

\begin{figure}[t]
\centering
   \includegraphics[width=0.35\linewidth,angle=-90,trim=1in 2in 1.5in 1in,clip]{Signconsistentwb.eps} \\
   \caption{Geometric interpretation of the conditions for sign consistency. Shaded area represents $\Omega_{A_0,\mathbf{s}_0}$,
   where $\mathbf{s}_0=\mbox{sgn}(\bm{\beta}_{0A_0})$.
    \label{fig:consistent}}
\end{figure}

Now we may recover some of the conditions for establishing consistency of the Lasso in the literature.
In what follows, let $\mathbf{W}=\mathbf{I}_p$ in \eqref{eq:muA0} and \eqref{eq:sigmaA0}. Condition (C2) is
the strong irrepresentable condition \citep{ZhaoYu06,Meinshausen06,Zou06}:
\begin{equation*}
\|\mathbf{C}_{I_0A_0}\mathbf{C}_{A_0A_0}^{-1}\mathbf{s}_{0}\|_{\infty} \leq c \in (0,1).
\end{equation*}
Assume that the minimum eigenvalue of $\mathbf{C}_{A_0A_0}$ is bounded from below by $\phi_0>0$, which is equivalent to condition (6) in \cite{ZhaoYu06}. Condition (C1) holds if
\begin{equation}\label{eq:C1}
\frac{\lambda\|\mathbf{C}_{A_0A_0}^{-1}\mathbf{s}_{0}\|_{\infty}}{\inf_{j \in A_0} |\beta_{0j}|}
\leq \frac{\lambda \phi_0^{-1} \|\mathbf{s}_0\|_2}{\inf_{j \in A_0} |\beta_{0j}|}
=\frac{\phi_0^{-1}\lambda \sqrt{q_0}}{\inf_{j \in A_0} |\beta_{0j}|}\to 0.
\end{equation}
This shows that some version of a beta-min condition is necessary to enforce
a lower bound for $\inf_{j \in A_0} |\beta_{0j}|$. For example, we may assume that
\begin{equation}\label{eq:betaminEA}
\lim_{n\to\infty} n^{a_1}\inf_{j \in A_0} |\beta_{0j}| \geq M_3,
\end{equation}
for some positive constants $M_3$ and $a_1$, which is the same as (8) in \cite{ZhaoYu06}.
Then one needs to choose $\lambda=o(n^{-a_1}/\sqrt{q_0})\to 0$ for (C1) to hold.

Let $\mathbf{d}=(d_1,\ldots,d_p)$ such that $\mathbf{d}_{A_0}=\sigma^2\mbox{diag}(\mathbf{C}_{A_0A_0}^{-1})$
and $\mathbf{d}_{I_0}=\sigma^2 \mbox{diag}(\mathbf{C}_{I_0\mid A_0})$.
To establish condition (C3), assume that $\bm{\varepsilon}\sim\mathcal{N}_n(\mathbf{0},\sigma^2 \mathbf{I}_n)$ and
$\mbox{diag}(\mathbf{C})$ is bounded from above.
Then all $d_j$ are bounded and let $d_*<\infty$ be an upper bound of $\{d_1,\ldots,d_p\}$.
Furthermore, $Z_j$ follows a univariate normal distribution:
For $j\in A_0$, $Z_j \sim \mathcal{N}(\mu^0_j,n^{-1}d_j)$ and for $j\in I_0$,
$Z_j \sim \mathcal{N}(\mu^0_j,n^{-1}\lambda^{-2} d_j)$ according to \eqref{eq:sigmaA0} with $\mathbf{W}=\mathbf{I}_p$.
By \eqref{eq:C1} and \eqref{eq:betaminEA},
$\delta_j=|\mu^0_j|\geq \frac{1}{2}M_3 n ^{-a_1}$ for $j\in A_0$ as $n\to \infty$, and
\begin{equation}\label{eq:C3A0}
P\left(\sup_{j \in A_0}|Z_j - \mu^0_j| \geq \frac{1}{2}M_3 n ^{-a_1} \right)
\leq 2 q_0 \exp\left(-\frac{M_3^2 n^{1-2a_1}}{8 d_*}  \right) \to 0,
\end{equation}
as long as $a_1<1/2$ and $q_0<n$. Since $\delta_j=1-c$ for all $j\in I_0$,
\begin{equation}\label{eq:C3I0}
P\left(\sup_{j \in I_0}|Z_j - \mu^0_j| \geq 1-c \right)
\leq  2 \exp\left(-\frac{n \lambda^2(1-c)^2}{2 d_*}+ \log p\right)\to 0,
\end{equation}
if $(\log p)/(n \lambda^2) \to 0$ and $n\lambda^2 \to \infty$. Clearly, the above two inequalities imply \eqref{eq:C3}.
Therefore, $\lambda$ must satisfy
$\sqrt{(\log p)/n} \ll \lambda = o(n^{-a_1}/\sqrt{q_0})$, which implies that $\sqrt{q_0 (\log p)/n} =o(n^{-a_1})$.
This is consistent with the beta-min condition in \cite{Meinshausen06}: $\inf_{j \in A_0} |\beta_{0j}|\gg \sqrt{q_0 (\log p)/n}$.
In summary, choosing $a_1,a_2,a_3>0$ such that $a_2+a_3< 1-2a_1$, the Lasso can be consistent for model selection
with $q_0 = O(n^{a_2})$ and $p=O(\exp(n^{a_3}))$, both diverging with $n$.
For more general scaling of $(n, p, q_0)$, see the work by \cite{Wainwright09}.

\begin{remark}
The term $\log p$ in \eqref{eq:C3I0} can be replaced by $\log(p-q_0)$, which will improve the bound
if $q_0/p$ does not vanish as $n\to\infty$. Moreover, both inequalities \eqref{eq:C3A0} and \eqref{eq:C3I0} are applicable
to sub-Gaussian noise.
\end{remark}

\subsection{Bayesian interpretation}\label{sec:Bayesian}

It is well-known that the Lasso can be interpreted as the mode of the posterior distribution of $\bm{\beta}$
under a Laplace prior. However, the posterior distribution itself is continuous on $\mathbb{R}^p$.
If we draw $\bm{\beta}$ from this posterior distribution, every component of $\bm{\beta}$ will be nonzero with
probability one. In this sense, sampling from this posterior distribution does not provide a direct solution
to model selection, which seems unsatisfactory from a Bayesian perspective. Here,
we discuss a different Bayesian interpretation of the Lasso-type estimator $\hat{\bmbeta}$
from a sampling distribution point of view.

Assume that $\text{rank}(\mathbf{X})=p< n$ and thus $\mathbf{C}$ is invertible.
Under the noninformative prior $p(\bm{\beta},\sigma^2) \propto 1/\sigma^2$ and
the assumption that $\bm{\varepsilon}\sim \mathcal{N}_n (\mathbf{0}, \sigma^2 \mathbf{I}_n)$, the conditional and marginal
posterior distributions of $\bm{\beta}$ are
\begin{eqnarray}
\bm{\beta} \mid \sigma^2, \bfy & \sim & \mathcal{N}_p(\hat{\bmbeta}^{\text{OLS}},n^{-1}\sigma^2\mathbf{C}^{-1}), \label{eq:condposterior}\\
\bm{\beta} \mid \bfy & \sim & t_{n-p}(\hat{\bmbeta}^{\text{OLS}},n^{-1}\hat{\sigma}^2\mathbf{C}^{-1}), \label{eq:marginalposterior}
\end{eqnarray}
where $\hat{\sigma}^2$ is given by \eqref{eq:estsigma2}
with $\check{\bmbeta}=\hat{\bmbeta}^{\text{OLS}}$ and $t_{n-p}(\bm{\mu},\bm{\Sigma})$ is the multivariate $t$ distribution
with $(n-p)$ degrees of freedom, location $\bm{\mu}$, and scale matrix $\bm{\Sigma}$.

Following the decision theory framework, let
$\bm{\eta}\in\mathbb{R}^p$ be a decision regarding $\bm{\beta}$ that incurs the loss
\begin{equation}\label{eq:Bayesloss}
\ell_{B}(\bm{\eta},\bm{\beta})=
\frac{1}{2}(\bm{\eta}-\bm{\beta})^{\mathsf{T}} \mathbf{C} (\bm{\eta}-\bm{\beta}) + \lambda \| \mathbf{W} \bm{\eta}\|_1.
\end{equation}
Since the covariance of $\bm{\beta}$ is proportional to $\mathbf{C}^{-1}$ with respect to
the posterior distribution \eqref{eq:condposterior} or \eqref{eq:marginalposterior},
$\ell_{B}(\bm{\eta},\bm{\beta})$ is essentially the squared Mahalanobis distance between $\bm{\eta}$
and $\bm{\beta}$, plus a weighted $\ell_1$ norm of $\bm{\eta}$ to encourage sparsity.
Denote by $\tilde{\bmbeta}$ the optimal decision that minimizes the
loss $\ell_B$ for a given $\bm{\beta}$, i.e.,
$\tilde{\bmbeta}=\arg\min_{\bm{\eta}} \ell_{B}(\bm{\eta},\bm{\beta})$.
Let $\tilde{\mathbf{S}}$ be the subgradient of $\|\bm{\eta}\|_1$ at $\tilde{\bmbeta}$.
The KKT condition for $\tilde{\bmbeta}$ is
\begin{equation}\label{eq:Bayesiangrad}
\mathbf{C} \tilde{\bmbeta} + \lambda \mathbf{W} \tilde{\mathbf{S}} =\mathbf{C}\bm{\beta}.
\end{equation}
Since $\bm{\beta}$ is a random vector in Bayesian inference, the distribution of $\bm{\beta}$
determines the joint distribution of $\tilde{\bmbeta}$ and $\tilde{\mathbf{S}}$ via the above KKT condition.
Represent $(\tilde{\bmbeta},\tilde{\mathbf{S}})$ by its equivalent form $(\tilde{\bmbeta}_{\tilde{\mathcal{A}}},\tilde{\mathbf{S}}_{\tilde{\mathcal{I}}},\tilde{\mathcal{A}})$ in the same way as for $(\hat{\bmbeta},\mathbf{S})$ in Section~\ref{sec:distribution}.

The conditional posterior distribution \eqref{eq:condposterior} implies that
$\mathbf{C}\bm{\beta} \mid \sigma^2, \bfy\sim \mathcal{N}_p(\mathbf{C}\hat{\bmbeta}^{\text{OLS}}, n^{-1}{\sigma^2}\mathbf{C})$.
Thus, conditional on $\bfy$ and $\sigma^2$, Equation \eqref{eq:Bayesiangrad} implies that
\begin{equation}\label{eq:Bayesiangraddistr}
\mathbf{C} \tilde{\bmbeta} + \lambda \mathbf{W} \tilde{\mathbf{S}} - \mathbf{C}\hat{\bmbeta}^{\text{OLS}} \eqinL \mathbf{U},
\end{equation}
where $\mathbf{U}\sim\mathcal{N}_p(\mathbf{0},n^{-1}{\sigma^2}\mathbf{C})$. One sees that \eqref{eq:Bayesiangraddistr}
is identical to the KKT condition \eqref{eq:lassoKKTinU} with $\hat{\bmbeta}^{\text{OLS}}$ in place of $\bm{\beta}$.
Therefore, the conditional distribution $[\tilde{\bmbeta}_{\tilde{\mathcal{A}}},\tilde{\mathbf{S}}_{\tilde{\mathcal{I}}},\tilde{\mathcal{A}}\mid \sigma^2, \bfy]$,
determined by \eqref{eq:Bayesiangraddistr},
is identical to the estimated sampling distribution $\hat{\pi}$ \eqref{eq:estjointdensity} under a normal error
distribution with $\bm{\beta}$ estimated by $\hat{\bmbeta}^{\text{OLS}}$, i.e., $\check{\bmbeta}=\hat{\bmbeta}^{\text{OLS}}$.
Furthermore, $\mathbf{C}\bm{\beta} \mid \bfy\sim t_{n-p}(\mathbf{C}\hat{\bmbeta}^{\text{OLS}}, n^{-1}\hat{\sigma}^2\mathbf{C})$ due to
\eqref{eq:marginalposterior}. By a similar reasoning, the conditional distribution
$[\tilde{\bmbeta}_{\tilde{\mathcal{A}}},\tilde{\mathbf{S}}_{\tilde{\mathcal{I}}},\tilde{\mathcal{A}}\mid \bfy]$ is the same as $\hat{\pi}$ if $\check{\bmbeta}=\hat{\bmbeta}^{\text{OLS}}$
and if $f_{\mathbf{U}}$ is estimated by the density of $t_{n-p}(\mathbf{0}, n^{-1}\hat{\sigma}^2\mathbf{C})$.
This motivates our proposal to use $t_{n-p}(\mathbf{0}, n^{-1}\sigma^2\mathbf{C})$
as a parametric model for ${\mathbf{U}}$ and estimate $\sigma^2$ from data
to construct $\hat{f}_{\mathbf{U}}$. The above discussion also provides a Bayesian justification for sampling from $\hat{\pi}$.

Under this framework, we may define
a point estimator $\hat{\bmbeta}^{\text{P}}=(\hat{\beta}^{\text{P}}_j)_{1:p}$ by the decision that
minimizes the posterior expectation of the loss $\ell_{B}(\bm{\eta},\bm{\beta})$,
\begin{equation}\label{eq:defBest}
\hat{\bmbeta}^{\text{P}}
\defi \arg\min_{\bm{\eta}}\int \ell_{B}(\bm{\eta},\bm{\beta}) p(\bm{\beta} \mid \bfy) d\bm{\beta},
\end{equation}
provided that the expectation exists. Although $\hat{\bmbeta}^{\text{P}}$ minimizes the posterior expected loss, its Bayes risk is not well-defined due to our use of an improper prior. To avoid any potential confusion, we call $\hat{\bmbeta}^{\text{P}}$ a posterior point estimator instead of a Bayes estimator. Taking subderivative of $\ell_{B}(\bm{\eta},\bm{\beta})$ with respect to $\bm{\eta}$ leads to the following equation to solve for the minimizer $\hat{\bmbeta}^{\text{P}}$:
\begin{equation}\label{eq:KKTBest}
\mathbf{C} \hat{\bmbeta}^{\text{P}} + \lambda \mathbf{W} \mathbf{S}^{\text{P}}
= \int  \mathbf{C}\bm{\beta} \cdot p(\bm{\beta} \mid \bfy) d\bm{\beta}=\mathbb{E}(\mathbf{C}\bm{\beta} \mid \bfy),
\end{equation}
where $\mathbf{S}^{\text{P}}$ is the subgradient of $\|\bm{\eta}\|_1$ at $\hat{\bmbeta}^{\text{P}}$. Under the noninformative prior,
the posterior mean $\mathbb{E}(\bm{\beta} \mid \bfy)=\hat{\bmbeta}^{\text{OLS}}$. In this case,
$\mathbb{E}(\mathbf{C}\bm{\beta} \mid \bfy)=n^{-1}\mathbf{X}^{\mathsf{T}}\bfy$ and Equation \eqref{eq:KKTBest}
is identical to the KKT condition \eqref{eq:lassograd} for the Lasso-type estimator $\hat{\bmbeta}$.
Therefore, $\hat{\bmbeta}$ can be interpreted as the estimator \eqref{eq:defBest} that minimizes the posterior expected loss.

\begin{remark}\label{rm:Bayes}
These results provide a Bayesian interpretation of the Lasso-type estimator $\hat{\bmbeta}$ and its sampling distribution.
Assume a normal error distribution with a given $\sigma^2$ and the noninformative prior.
The posterior distribution of the optimal decision, $[\tilde{\bmbeta}\mid \bfy]$, is identical to the sampling distribution of $\hat{\bmbeta}$
assuming $\hat{\bmbeta}^{\text{OLS}}$ is the true coefficient vector. Therefore, a posterior probability interval for $\tilde{\bmbeta}$, the optimal decision, constructed according to $[\tilde{\bmbeta}\mid \bfy]$ is the same as the confidence interval
constructed according to $\hat{\pi}$ with
$\check{\bmbeta}=\hat{\bmbeta}^{\text{OLS}}$. Point estimation about $\bm{\beta}$ also coincides between the Bayesian
and the penalized least-squares methods ($\hat{\bmbeta}^{\text{P}}=\hat{\bmbeta}$).
Lastly, if we set $\lambda=0$ in the loss \eqref{eq:Bayesloss}, then the optimal decision $\tilde{\bmbeta}$ is simply $\bm{\beta}$.
In this special case, the aforementioned coincidences become the familiar correspondence between
the posterior distribution \eqref{eq:condposterior} and the sampling distribution of $\hat{\bmbeta}^{\text{OLS}}$
and that between the posterior mean and $\hat{\bmbeta}^{\text{OLS}}$.
\end{remark}

It is worth mentioning that, in a loose sense, this Bayesian interpretation also applies when $p>n$.
In this case, the posterior distribution \eqref{eq:condposterior} does not exist, but $[\mathbf{C}\bm{\beta} \mid\sigma^2, \bfy]$
is a well-defined normal distribution in $\text{row}(\mathbf{X})$. From KKT conditions \eqref{eq:Bayesiangrad} and \eqref{eq:KKTBest},
we see that the posterior point estimator $\hat{\bmbeta}^{\text{P}}$ and the posterior distribution $[\tilde{\bmbeta}\mid \bfy]$
only depend on $\mathbf{C}\bm{\beta}$. Therefore, they are well-defined and have the same coincidence
with the Lasso-type estimator and its sampling distribution.

\subsection{Bootstrap versus Monte Carlo}\label{sec:limitation}

We have demonstrated that Monte Carlo sampling via estimator augmentation has substantial advantages in approximating tail probabilities and conditional distributions, say $[\hat{\bmbeta}_A \mid \mathcal{A}=A]$, over direct sampling (or  bootstrap). The MH Lasso sampler also showed some improvement in efficiency when compared against direct sampling in the low-dimensional setting. Now we discuss some limitations of estimator augmentation relative to bootstrap.

The joint density of the augmented estimator is derived for a given $\lambda$, and thus does not take into account the randomness in $\lambda$ when it is chosen via a data-dependent way, say via cross-validation. Denote by $\hat{\bmbeta}(\bfy,\hat{\lambda}(\bfy))$ the Lasso-type estimator when $\lambda=\hat{\lambda}(\bfy)$, where $\hat{\lambda}(\bfy)$ is estimated from the data $\bfy$. We stress that the density in Theorem~\ref{thm:joint} or Theorem~\ref{thm:jointhigh} does not apply to the sampling distribution of $\hat{\bmbeta}(\bfy,\hat{\lambda}(\bfy))$ and it is only valid for $\hat{\bmbeta}(\bfy,\lambda)$ with $\lambda$ being fixed during the repeated sampling of $\bfy$. However, the direct sampler (or bootstrap in a similar way) can handle data-dependent $\lambda$ by adding one additional step to determine $\hat{\lambda}(\bfy^{(t)})$ after each draw of $\bfy^{(t)}$ in Routine~\ref{rt:directsample}.

Bootstrap and the direct sampler can be parallelized. The importance sampling algorithm (Routine~\ref{rt:IS}) can easily be parallelized as well, since it uses the direct sampler to generate proposals and calculates importance weights independently for each sample. An MCMC algorithm needs a certain number of burn-in iterations before the Markov chain reaches its stationary distribution. It seems that naively running multiple short chains in parallel may impair the overall efficiency due to the computational waist of multiple burn-in iterations. Initialized with one draw from the direct sampler, a Markov chain simulated by Routine~\ref{rt:dirmcmc}, however, reaches its equilibrium at the first iteration and thus is suitable for parallel computing. Its efficiency relative to direct sampling when both are parallelized can be calculated as follows.

Suppose our goal is to estimate $\mathbb{E}_{\pi}[g(\hat{\bmbeta})]$ and assume that $\text{Var}_{\pi}[g(\hat{\bmbeta})]=1$ without loss of generality. Assume that the time to run one iteration of the direct sampler allows for running $m$ iterations of an MCMC algorithm. Suppose that we have access to $K$ computing nodes and the available computing time from each node allows for the simulation of $(1+N_1)$ samples from the direct sampler, where $N_1$ may be small. Thus, on a single node we can run $N_2=mN_1$ MCMC iterations plus an initial draw from the direct sampler in the same amount of time. In other words, we can run Routine~\ref{rt:dirmcmc} for $1+N_2$ iterations to draw $\bm{\beta}^{(t)}$ for $t=1,\ldots,1+N_2$. Note that this Markov chain reaches equilibrium from $t=1$. Let $\rho_t=\text{cor}(g(\bm{\beta}^{(1)}),g(\bm{\beta}^{(t+1)}))$ and
\begin{equation*}
\psi(N)=1+2 \sum_{t=1}^{N-1} \left(1-\frac{t}{N}\right)\rho_t
\end{equation*}
for an integer $N\geq 1$. Then we have
\begin{equation*}
\text{Var}\left[\frac{1}{N_2+1} \sum_{t=1}^{N_2+1}g(\bm{\beta}^{(t)})\right]=\frac{1}{N_2+1}\psi(N_2+1)\defi V_2(N_2+1).
\end{equation*}
Denote by $V_1(N)=1/N$ the variance in estimating $g$ by the mean of an i.i.d. sample of size $N$, and let
\begin{equation*}
\gamma=\lim_{N\to\infty}\frac{V_1(N)}{V_2(mN)}= \frac{m}{\psi(\infty)}.
\end{equation*}
The efficiency of Routine~\ref{rt:dirmcmc} relative to direct sampling is
\begin{align*}
\frac{V_1(N_1+1)}{V_2(N_2+1)} &=\frac{mN_1+1}{N_1+1}\frac{1}{\psi(N_2+1)} \\
& >\frac{N_1}{N_1+1}\frac{m}{\psi(N_2+1)} \geq \frac{N_1}{N_1+1}\gamma,
\end{align*}
where we have assumed that $\psi(N_2+1)\leq \psi(\infty)$ for the last inequality. This assumption holds if $\psi(N)$ is nondecreasing in $N$. This derivation shows that Routine~\ref{rt:dirmcmc} will be more efficient than the direct sampler on each computing node if $N_1\geq 1/(\gamma-1)$, which can be as small as 1 when $\gamma>2$.
We have observed two decay patterns of the autocorrelation $\rho_t$ of the MLS in the simulation study in Section~\ref{sec:numerical}. For some components of $\hat{\bmbeta}$, $\rho_t$ is always positive before it decays to zero, in which case $\psi(N)$ is obviously nondecreasing. For other components, $\rho_t$ first decreases monotonely to zero and then shows small fluctuations around zero. In the second case, we empirically observed that $\psi(N)$ is nondecreasing as well. The efficiency comparison in Table~\ref{tab:MSEjoint}, with $N_1$ and $N_2=mN_1$ both large, suggests that for most functions estimated there, $\gamma\in(2,3)$  for datasets A and B and $\gamma\in (1.2,1.6)$ for the other two datasets. Therefore, as long as we need to run a few iterations of the direct sampler on each node, parallelizing Routine~\ref{rt:dirmcmc} can bring computational gain. Of course, if the number of computing nodes $K$ is so large that only one draw is needed from each node, direct sampling or bootstrap will be a better choice.

\subsection{Concluding remarks}\label{sec:discuss}

Utilizing the density of an augmented estimator, this article develops
MCMC and IS methods to approximate sampling distributions
in $\ell_1$-penalized linear regression. This approach is clearly different from
existing methods based on resampling or asymptotic approximation. The numerical results have already
demonstrated the substantial gain in efficiency and the great flexibility offered by this approach.
These results are mostly for a proof of principle,
and there is room for further development of more efficient Monte Carlo algorithms based on
the densities derived in this article.

In principle, the idea of estimator augmentation can be applied to the use of concave penalties in linear regression
\citep{FrankFriedman93,FanLi01,Friedman08,Zhang10} for studying the sampling distribution.
However, there are at least two additional technical
difficulties for the high-dimensional setting. First, we need to find conditions for the uniqueness of a concave-penalized
estimator in order to construct a bijection between $\mathbf{U}$ and the augmented estimator. Second, the constraint in
\eqref{eq:spacehigh} will become nonlinear in general, even for a fixed $\mathbf{s}_A$, when a concave penalty is used,
which means that the sample space is composed of a finite number of manifolds.
Another future direction is to investigate theoretically and empirically the finite-sample
performance in variable selection by the Lasso sampler which may take into account the uncertainty in parameter estimation in a coherent way.

\section{Proofs}\label{sec:proofs}

Let $n\geq 1$ and $p\geq 2$ throughout this section.

\subsection{Proof of Theorem~\ref{thm:threshold}}

\begin{lemma}\label{lm:truncated}
Let $\mathbf{Z}\in\mathbb{R}^p$ be a random vector, $K \in \mathscr{R}^p$, and $\mathbf{Z}_K$ be the truncation of $\mathbf{Z}$ to $K$ such that $P(\mathbf{Z}_K\in B)=P(\mathbf{Z}\in B \mid \mathbf{Z} \in K)$ for $B \in \mathscr{R}^p$. If $P(\mathbf{Z} \in K)\geq 1-\alpha>0$, then
\begin{align*}
\sup_{B\in \mathscr{R}^p}|P(\mathbf{Z}_K\in B)-P(\mathbf{Z} \in B)|\leq \alpha.
\end{align*}
\end{lemma}
\begin{proof}
For any $B\in\mathscr{R}^p$, $P(\mathbf{Z} \in B \cap K) = P(\mathbf{Z}_K\in B) P(\mathbf{Z}\in K)$ and thus
\begin{align*}
0 \leq P(\mathbf{Z}_K\in B)- P(\mathbf{Z} \in B \cap K) = P(\mathbf{Z}_K\in B) P(\mathbf{Z}\in K^c) \leq \alpha.
\end{align*}
On the other hand,
\begin{align*}
0 \leq P(\mathbf{Z}\in B)- P(\mathbf{Z} \in B \cap K) =P(\mathbf{Z}\in B \cap K^c) \leq \alpha.
\end{align*}
Therefore, $|P(\mathbf{Z}_K\in B)-P(\mathbf{Z} \in B)|\leq \alpha$ for any $B$ and the conclusion follows.
\end{proof}

\begin{lemma}\label{lm:Veq}
Assume that $\check{\bmbeta}$ satisfies conditions (i) and (ii) in Lemma~\ref{lm:fixknownvar}, and let \eqref{eq:assbeta0j} be satisfied.
Then $V(\bm{\delta}; \check{\bmbeta}, \mathbf{u})=V(\bm{\delta}; \bm{\beta}_0, \mathbf{u})$ for any $\mathbf{u} \in \mathbb{R}^p$ if $\|\bm{\delta}\|_{\infty}\leq M_1$.
\end{lemma}
\begin{proof}
The assumptions on $\check{\bmbeta}$ imply that $\mathcal{A}(\check{\bmbeta})=A_0$ and $\mbox{sgn}(\check{\beta}_j)=\mbox{sgn}(\beta_{0j})$ for all $j \in A_0$. By the definition of $V$ \eqref{eq:defVn} it then suffices to show that
\begin{align}\label{eq:signequal}
| \beta_{0j}+ r_n^{-1}\delta_{j}|- | \beta_{0j}|=| \check{\beta}_{j}+ r_n^{-1}\delta_{j}|- | \check{\beta}_{j}|
\end{align}
for all $j\in A_0$. By the definition of $\eta$ in \eqref{eq:difftdbeta}, for $j\in A_0$
\begin{align*}
\eta |\beta_{0j}| \geq |\check{\beta}_{j}-\beta_{0j}| \geq |\beta_{0j}| - |\check{\beta}_{j}|
\end{align*}
and therefore
\begin{align*}
|\check{\beta}_{j}| \geq (1-\eta) |\beta_{0j}| > M_1/r_n \geq |r_n^{-1}\delta_j|,
\end{align*}
where we have used \eqref{eq:assbeta0j} and that $\|\bm{\delta}\|_{\infty}\leq M_1$. Consequently, for $j \in A_0$ we have
\begin{align*}
| \check{\beta}_{j}+ r_n^{-1}\delta_{j}|- | \check{\beta}_{j}|= \mbox{sgn}(\check{\beta}_j) \delta_j/r_n.
\end{align*}
On the other hand, by \eqref{eq:assbeta0j} and $\eta \in[0,1)$,
$|\beta_{0j}| > M_1/r_n \geq |r_n^{-1}\delta_j|$ for $j\in A_0$ and thus
\begin{align*}
| \beta_{0j}+ r_n^{-1}\delta_{j}|- | \beta_{0j}|= \mbox{sgn}(\beta_{0j}) \delta_j/r_n.
\end{align*}
Now \eqref{eq:signequal} follows since $\mbox{sgn}(\check{\beta}_j)=\mbox{sgn}(\beta_{0j})$ for all $j \in A_0$.
\end{proof}

\begin{proof}[Proof of Lemma~\ref{lm:fixknownvar}]
Define $\tilde{\bmdelta}=\arg\min_{\bm{\delta}}V(\bm{\delta};\bm{\beta}_0,\mathbf{U}^*)$, which follows the same distribution as $\hat{\bm{\delta}}$ \eqref{eq:defdelta}. This is because $\mathbf{U}^*\eqinL \mathbf{U}$ by fixing $\hat{\sigma}=\sigma$ and $V(\bm{\delta};\bm{\beta}_0,\mathbf{u})$ has a unique minimizer for any $\mathbf{u}$ if the columns of $\mathbf{X}$ are in general position (Lemma~\ref{lm:prelim}). Consequently,
\begin{align*}
\nu[\hat{\bm{\delta}} \mid \|\hat{\bm{\delta}}\|_{\infty} < M_1] =\nu[\tilde{\bmdelta}\mid \|\tilde{\bmdelta}\|_{\infty} < M_1].
\end{align*}
Let $\mathcal{K}=\{\bm{\delta}\in\mathbb{R}^p:\|\bm{\delta}\|_{\infty}< M_1\}$.
According to Lemma~\ref{lm:Veq}, $V(\bm{\delta};\check{\bmbeta},\mathbf{U}^*)=V(\bm{\delta};\bm{\beta}_0,\mathbf{U}^*)$ for all $\bm{\delta}\in\mathcal{K}$.
As the unique minimizer of $V(\bm{\delta};\check{\bmbeta},\mathbf{U}^*)$ \eqref{eq:defdeltastar}, $\|\bmdelta^* \|_{\infty}< M_1$ implies that $\bmdelta^*$ is also a local minimizer of $V(\bm{\delta};\bm{\beta}_0,\mathbf{U}^*)$. Since $V(\bm{\delta};\bm{\beta}_0,\mathbf{U}^*)$ is convex in $\bm{\delta}$ and has only a unique minimizer $\tilde{\bmdelta}$, we must have $\bmdelta^*=\tilde{\bmdelta}$ and $\|\tilde{\bmdelta} \|_{\infty}< M_1$. Furthermore, using the same argument in the other direction, one can show that $\|\tilde{\bmdelta} \|_{\infty}< M_1$ implies $\|\bmdelta^* \|_{\infty}< M_1$, and thus $\{\|\tilde{\bmdelta} \|_{\infty}< M_1\}$ is equivalent to $\{\|\bmdelta^* \|_{\infty}< M_1\}$. This completes the proof.
\end{proof}

\begin{proof}[Proof of Theorem~\ref{thm:threshold}]
Let $E_1$ be the event that $\check{\bmbeta}$ satisfies conditions (i) and (ii) in Lemma~\ref{lm:fixknownvar} and $E_2=\{\|\hat{\bm{\delta}}\|_{\infty} < M_1\}$. We first show that \eqref{eq:supPrdiff} holds on $E_1$. Obviously, \eqref{eq:assbeta0j} holds because of \eqref{eq:betamin}. The argument in the proof of Lemma~\ref{lm:fixknownvar} implies that, on event $E_1$,
\begin{align}
P(\|\bmdelta^*\|_{\infty} < M_1\mid \check{\bmbeta}) &=P(\|\tilde{\bmdelta}\|_{\infty} < M_1) \nonumber \\
& = P(E_2)\geq 1-\alpha_1, \label{eq:truncatePr}
\end{align}
where the second equality is due to $\tilde{\bmdelta}\eqinL \hat{\bm{\delta}}$. Let $\bm{\delta}^*_{\mathcal{K}}$ and $\hat{\bm{\delta}}_{\mathcal{K}}$ be the respective truncations of $\bmdelta^*$ and $\hat{\bm{\delta}}$ to $\mathcal{K}$. Lemma~\ref{lm:fixknownvar} implies $\nu[\bm{\delta}^*_{\mathcal{K}} \mid \check{\bmbeta}]=\nu[\hat{\bm{\delta}}_{\mathcal{K}}]$ on event $E_1$. A direct consequence is that on $E_1$, $P(\bm{\delta}^*_{\mathcal{K}}\in B \mid \check{\bmbeta})=P(\hat{\bm{\delta}}_{\mathcal{K}} \in B)$ for any $B\in\mathscr{R}^p$ and therefore,
\begin{align*}
\sup_{B\in \mathscr{R}^p} |P(\bmdelta^* \in B \mid \check{\bmbeta}) - P(\hat{\bm{\delta}} \in B)| \leq 2\alpha_1
\end{align*}
by Lemma~\ref{lm:truncated} and \eqref{eq:truncatePr}.

Next we find a lower bound for $P(E_1)$. Since $\inf_{A_0}|\beta_{0j}|> 2\tau$ \eqref{eq:betamin}, on event $\mathcal{E}$, we have $\mathcal{A}(\check{\bmbeta})=A_0$, according to Lemma~\ref{lm:screen}. By construction $\check{\bmbeta}_{A_0}=\hat{\bmbeta}_{A_0}$ \eqref{eq:thresholddef} and consequently
\begin{equation}\label{eq:etaupper}
\eta =  \sup_{j \in A_0} \frac{|\hat{\beta}_{j}-\beta_{0j}|}{| \beta_{0j}|} \leq \frac{\|\hat{\bmbeta}-\bm{\beta}_{0}\|_2}{\inf_{A_0}| \beta_{0j}|} < \frac{1}{2}
\end{equation}
again on $\mathcal{E}$. Thus, $P(E_1)\geq P(\mathcal{E}) \geq 1-2p^{1-u}$.
\end{proof}

\subsection{Proof of Theorem~\ref{thm:thresholdsigmaest}}

\begin{lemma}\label{lm:devVmin}
Let $\bm{\gamma}$ be any minimizer of $V(\bm{\delta};\check{\bmbeta},\mathbf{u})$ for $\check{\bmbeta}\in\mathbb{R}^p$ and $\mathbf{u} \in \mathbb{R}^p$. For any $\bm{\Delta} \in \mathbb{R}^p$, we have
\begin{align}\label{eq:bounddevVmin}
V(\bm{\gamma}+\bm{\Delta};\check{\bmbeta},\mathbf{u}) \geq V(\bm{\gamma};\check{\bmbeta},\mathbf{u}) + \frac{n}{2r_n^2} \bm{\Delta}^{\mathsf{T}} \mathbf{C} \bm{\Delta}.
\end{align}
\end{lemma}
\begin{proof}
Let $\mathbf{b}=r_n^{-1}\bm{\gamma}+\check{\bmbeta}=(b_j)_{1:p}$ and $\bm{\Delta}=(\Delta_j)_{1:p}$. Direct calculations give
\begin{align*}
& \quad V(\bm{\gamma}+\bm{\Delta};\check{\bmbeta},\mathbf{u})- V(\bm{\gamma};\check{\bmbeta},\mathbf{u})  \\
 & = \frac{n}{2r_n^2} \bm{\Delta}^{\mathsf{T}} \mathbf{C} (\bm{\Delta} + 2 \bm{\gamma})
- \frac{n}{r_n}\bm{\Delta}^{\mathsf{T}} \mathbf{u}
+ n\lambda \sum_{j=1}^p w_j (| b_j+ r_n^{-1}\Delta_{j}|-|b_j|).
\end{align*}
The KKT condition for $\bm{\gamma}$ to minimize $V(\bm{\delta};\check{\bmbeta},\mathbf{u})$ is
\begin{align}\label{eq:KKTminV}
\mathbf{C}(\mathbf{b}-\check{\bmbeta}) + \lambda \mathbf{W} \mathbf{s} -\mathbf{u} =\mathbf{0},
\end{align}
where $\mathbf{s}=(s_j)_{1:p}$ is the subgradient of $\|\bm{\beta}\|_1$ at $\mathbf{b}$. By the definition of a subgradient,\begin{align}\label{eq:subgradineq}
| b_j+ r_n^{-1}\Delta_{j}|-|b_j| \geq s_j r_n^{-1}\Delta_j
\end{align}
for all $j=1,\ldots,p$. Now we have
\begin{align*}
& V(\bm{\gamma}+\bm{\Delta};\check{\bmbeta},\mathbf{u})-V(\bm{\gamma};\check{\bmbeta},\mathbf{u}) \\
& \quad \geq \frac{n}{2r_n^2} \bm{\Delta}^{\mathsf{T}} \mathbf{C} \bm{\Delta} + \frac{n}{r_n^2} \bm{\Delta}^{\mathsf{T}} \mathbf{C} \bm{\gamma}
- \frac{n}{r_n}\bm{\Delta}^{\mathsf{T}} \mathbf{u} +\frac{n\lambda}{r_n} \sum_{j=1}^p w_j s_j \Delta_j \\
& \quad = \frac{n}{2r_n^2} \bm{\Delta}^{\mathsf{T}} \mathbf{C} \bm{\Delta} + \frac{n}{r_n} \bm{\Delta}^{\mathsf{T}}
\left[\mathbf{C}(\mathbf{b}-\check{\bmbeta}) - \mathbf{u} + \lambda \mathbf{W} \mathbf{s} \right] = \frac{n}{2r_n^2} \bm{\Delta}^{\mathsf{T}} \mathbf{C} \bm{\Delta},
\end{align*}
where we have used \eqref{eq:subgradineq} and \eqref{eq:KKTminV}.
\end{proof}

\begin{lemma}\label{lm:scaleu}
Assume that $|\mathcal{A}(\check{\bmbeta})|\leq q$. Let $c>0$, $\mathbf{u}=(u_j)_{1:p}\in\mathbb{R}^p$, and $\bm{\gamma}$ be any minimizer of $V(\bm{\delta};\check{\bmbeta},\mathbf{u})$. If Assumption RE$(q,3)$ is satisfied and $|u_j|\leq w_j\lambda/2$ for all $j=1,\ldots,p$, then
\begin{align}\label{eq:boundscaleu}
|V(\bm{\gamma};\check{\bmbeta},c\mathbf{u})-V(\bm{\gamma};\check{\bmbeta},\mathbf{u})|
\leq \frac{8 w_{\max}|1-c|}{w_{\min}\kappa^2(q,3)} n\lambda^2 \sum_{A_0} w_j^2.
\end{align}
\end{lemma}
\begin{proof}
Let $\mathbf{b}=r_n^{-1}\bm{\gamma}+\check{\bmbeta}$. Direct calculations give
\begin{align*}
 |V(\bm{\gamma};\check{\bmbeta},c\mathbf{u})-V(\bm{\gamma};\check{\bmbeta},\mathbf{u})|
& =  n \left|(1-c) \mathbf{u}^{\mathsf{T}} (\mathbf{b}-\check{\bmbeta}) \right| \\
& \leq  n |1-c| \cdot\|\mathbf{u}\|_{\infty} \|\mathbf{b}-\check{\bmbeta} \|_1.
\end{align*}
It is seen from \eqref{eq:KKTminV} that $\mathbf{b}$ is a minimizer of the loss \eqref{eq:lassoloss} if $\check{\bmbeta}$ is the true coefficient vector and if $\mathbf{X}^{\mathsf{T}}\bm{\varepsilon}/n=\mathbf{u}$. Inequality \eqref{eq:boundL1} in Lemma~\ref{lm:lassobounds} applied under these assumptions leads to
\begin{align*}
\|\mathbf{b}-\check{\bmbeta} \|_1 \leq \frac{16\lambda}{\kappa^2(q,3)} \sum_{A_0} \frac{w_j^2}{w_{\min}}
\end{align*}
if $|u_j|\leq w_j\lambda/2$ for all $j$. Moreover, $\|\mathbf{u}\|_{\infty} \leq w_{\max} \lambda/2$ and
hence \eqref{eq:boundscaleu} follows.
\end{proof}

\begin{proof}[Proof of Lemma~\ref{lm:fixunknownvar}]
To simplify notation, let $\hat{c}=\sigma/\hat{\sigma}$ and
\begin{equation*}
h=\frac{8 w_{\max}\zeta}{w_{\min}\kappa^2(q,3)} n\lambda^2 \sum_{A_0} w_j^2.
\end{equation*}
If $|U_j^*| \leq (1-\zeta)w_j\lambda/2$, then $|U_j^*| \leq w_j\lambda/2$ and by \eqref{eq:diffhsigma} $|\hat{c}U_j^*| \leq w_j\lambda/2$. Lemma~\ref{lm:scaleu} with \eqref{eq:diffhsigma} implies
\begin{align*}
|V(\bmdelta^*;\check{\bmbeta},\hat{c}\mathbf{U}^*)-V(\bmdelta^*;\check{\bmbeta},\mathbf{U}^*)| &\leq h,  \\
|V(\bmdelta^*_0;\check{\bmbeta},\mathbf{U}^*)-V(\bmdelta^*_0;\check{\bmbeta},\hat{c}\mathbf{U}^*)| &\leq h.
\end{align*}
Let $\bm{\Delta}=\bmdelta^*-\bmdelta^*_0$. Now we have
\begin{align*}
V(\bmdelta^*;\check{\bmbeta},\mathbf{U}^*) & \geq V(\bmdelta^*;\check{\bmbeta},\hat{c}\mathbf{U}^*) - h \\
& \geq V(\bmdelta^*_0;\check{\bmbeta},\hat{c}\mathbf{U}^*) + \frac{n}{2r_n^2} \bm{\Delta}^{\mathsf{T}} \mathbf{C} \bm{\Delta} -h \\
& \geq V(\bmdelta^*_0;\check{\bmbeta},\mathbf{U}^*)+ \frac{n}{2r_n^2} \bm{\Delta}^{\mathsf{T}} \mathbf{C} \bm{\Delta} -2h,
\end{align*}
where the second inequality is due to Lemma~\ref{lm:devVmin}. Lastly, since $V(\bmdelta^*_0;\check{\bmbeta},\mathbf{U}^*)\geq V(\bmdelta^*;\check{\bmbeta},\mathbf{U}^*)$ by definition \eqref{eq:defdeltastar}, $\bm{\Delta}^{\mathsf{T}} \mathbf{C} \bm{\Delta} \leq 4 r_n^2 h/n$ which coincides with \eqref{eq:boundquadratic}.
\end{proof}

\begin{proof}[Proof of Theorem~\ref{thm:thresholdsigmaest}]
Recall that $\mathcal{E}$ is the event $\cap_{j=1}^p\{|U_j| \leq w_j\lambda/2\}$. Since the distribution of $\bmdelta^*_0$ does not depend on $\hat{\sigma}$ and is identical to the distribution of $\bmdelta^*$ when $\hat{\sigma}$ is fixed to the true noise level $\sigma$, \eqref{eq:supPrdiffdelta0} follows immediately from \eqref{eq:supPrdiff} which holds on $\mathcal{E}$.

Let $E_3$ be the event in \eqref{eq:diffhsigma} and $\mathcal{E}^*$ be the event that $\cap_{j=1}^p\{|U_j^*| \leq (1-\zeta)w_j\lambda/2\}$.
By Lemma~\ref{lm:screen}, on $\mathcal{E}$ we have $\mathcal{A}(\check{\bmbeta})=A_0$, $|\mathcal{A}(\check{\bmbeta})|=q_0\leq q$ and by \eqref{eq:etaupper}
\begin{equation}\label{eq:bmbetacheck}
\inf_{A_0} |\check{\beta}_j| \geq \frac{1}{2} \inf_{A_0} |\beta_{0j}| > \tau.
\end{equation}
Therefore, all the assumptions on $\check{\bmbeta}$ and $\hat{\sigma}$ in Lemma~\ref{lm:fixunknownvar} are satisfied on $\mathcal{E} \cap E_3$, which happens with probability at least $1-(\alpha_2+2p^{1-u})$. Moreover, the conditional probability of \eqref{eq:boundquadratic} given $(\check{\bmbeta},\hat{\sigma})$ is at least
\begin{equation*}
P(\mathcal{E}^* \mid\hat{\sigma}) \geq 1-2p^{1-u},
\end{equation*}
by choosing $\lambda \geq (\hat{\sigma}/\sigma)\lambda_0/(1-\zeta)$. For the lower bound of the above probability, see \eqref{eq:condlambda} in Lemma~\ref{lm:lassobounds} with $(1-\zeta)w_j$ in place of $w_j$ and $\hat{\sigma}$ in place of $\sigma$. As $\hat{\sigma}/\sigma \leq 1+\zeta$ on $E_3$, it suffices to choose $\lambda$ as in \eqref{eq:condlambdasigmaest}. What remains is to show that $\bm{\Delta}=\bmdelta^*-\bmdelta^*_0$ is $M_2$-sparse on the event $\mathcal{E}^*$. Then \eqref{eq:L2Deltabound} follows from \eqref{eq:boundquadratic} and the definition of $\phi_{\min}(M_2)$ \eqref{eq:remin}.  Regarding $\check{\bmbeta}$ and $\hat{\sigma}$ as the true parameters, Lemma~\ref{lm:screen} with \eqref{eq:bmbetacheck} implies that $\mathcal{A}(\check{\bmbeta}) \subseteq \mathcal{A}({\bmbeta}^*)$ on $\mathcal{E}^*$ and therefore, $|\mathcal{A}(\bmdelta^*)| \leq |\mathcal{A}({\bmbeta}^*)|\leq M_2/2$ by \eqref{eq:boundsize}. Since $\mathcal{E}^*$ with \eqref{eq:diffhsigma} implies $|(\sigma/\hat{\sigma})U_j^*| \leq w_j\lambda/2$ for all $j$, by a similar reasoning we also have $|\mathcal{A}(\bmdelta^*_0)| \leq M_2/2$ and thus $|\mathcal{A}(\bm{\Delta})|\leq M_2$ on $\mathcal{E}^*$.
\end{proof}

redpajama_set_name: RedPajamaArXiv