\section{Introduction} \label{sec: intro}
The study of random networks has been a topic of great interest in recent years, e.g. see~\cite{kolaczyk2009statistical} and~\cite{newman2010networks}. A network is defined as a structure composed of nodes and edges connecting nodes in various relationships.
The observed network can be represented by an $N\times N$ adjacency matrix $\bm{Y} = (Y_{ij})_{i,j = 1,\ldots,N}$, where $N$ is the total number of nodes within the network. For a binary relation network, as considered here, $Y_{ij}=1$ if there is an edge from node $i$ to node $j,$ and $0$ otherwise. In the following we identify an adjacency matrix $\bm{Y}$ with the network itself.

Most relational phenomena are dependent phenomena, and dependence is often of substantive interest. \cite{frank1986markov} and  \cite{wasserman1996logit} introduced exponential random graph models which allow the modelling of a wide range of dependences of substantive interest, including transitive closure. For such models, $Y_{ij}\in\{0,1\}$ and the distribution of $\bm{Y}$ is assumed to follow the exponential family form
$P_{\bm{\theta}}(\bm{Y}=\bm{y})=\exp\left(\bm{\theta}\cdot \bm{T}(\bm{y})-\phi(\bm{\theta})\right), \bm{y}\in\mathcal{Y},$
where
$\phi(\bm{\theta})=-\log\left(\sum_{\bm{y}\in\mathcal{Y}}\exp(\bm{\theta}\cdot \bm{T}(\bm{y}))\right)$ and $\bm{T}(\bm{y}):\mathcal{Y}\rightarrow \mathbb{R}^q,$ are the sufficient statistics, e. g. the total number of edges.
However, as mentioned in \cite{schweinberger2014local}, exponential random graph models are lacking neighborhood structure, and that makes modelling dependencies challenging for such networks. Neighborhoods (or communities, blocks) are in general defined as a group of individuals (nodes), such that individuals within a group interact with each other more frequently than with those outside the group. Recently, \cite{schweinberger2014local} proposed the concept of local dependence in stochastic networks. This concept allows for dependence within neighborhoods, while different neighborhoods are independent.

In contrast to that, our work is considering dependence between blocks, while the connections within blocks are assumed independent. We also assume the blocks to be known. We then propose to analyze dependencies between blocks by means of graphical models. To this end, we assume an undirected network so that
\begin{align} \label{model: level 2}
	Y_{ij}|(\bm{P},\bm{z})\sim {\rm Bernoulli}\left(p_{\bm{z}[i],\bm{z}[j]}\right),
\end{align}
where $\bm{z}[i] \in \wb{K} := \{1, \cdots, K\}, i = 1,\ldots, N$ indicate block memberships in one of $K$ blocks; $p_{k,\ell}, k,\ell \in \wb{K}$ govern the intensities of the connectivities within and between blocks, $0 < p_{k,l} < 1$; and $\bm{P} = (p_{k, l})_{k, l\in\wb{K}}$ is a $K\times K$ symmetric matrix. We then put a Gaussian logistic model on the $p_{k,\ell}$. More precisely, for the diagonal elements $(p_{k,k})_{1\leq k\leq K}$, assume that
\begin{align} \label{model: level 1}
	\log\left(\frac{p_{k,k}}{1-p_{k,k}}\right)=\bm{x}_{k}^T\bm{\beta}+\epsilon_{k},\; 1\leq k \leq K,
\end{align}
where $\bm{x}_{k}$ is a $(L\times1)$ vector of given co-variables corresponding to block $k$, and $\bm{\beta}$ is the $(L\times1)$ parameter vector. Furthermore, $\bm{\epsilon} = (\epsilon_1,\cdots,\epsilon_K)^T$ with
\begin{align} \label{model: level1_error}
	\bm{\epsilon}\sim N(\bm{0},\bm{\Sigma}),
\end{align}
where $\bm{\Sigma}=(\sigma_{kl})_{1\leq k,l\leq K}$ is an nonsingular covariance matrix.
Each off-diagonal element $p_{k,l}$ ($k\neq l$) is assumed to be independent with all the other elements of $\bm{P}$. The latter assumption is made to simplify the exposition. A similar model can be found in \cite{xu2014dynamic}.

The dependence between the $p_{k,k}$ induces dependence between blocks. We can thus analyze this induced dependence in our network model, by using methods from Gaussian graphical models, via selecting the zeros in the precision matrix $\bm{\Sigma}^{-1}$. Adding dependencies between the $p_{k,\ell}$ with $k \ne \ell$ would increase the dimension of $\bm{\Sigma}$, and induce `second order dependencies' to the network structure, namely, dependencies of block connections between different pairs of blocks.

It is crucial to observe that this Gaussian graphical model is defined in terms of the $p_{k,k}$ (or, more precisely, in terms of their log-adds ratios), and that these quantities obviously are not observed. Thus, they need to be estimated from our network data, and, to this end, we here assume the availability of iid observations of the network. This estimation, in turn, induces {\em additional randomness} to our analysis of the graphical model. We are therefore facing similar challenges as in the analysis of Gaussian graphical models {\em under uncertainty}. However, our situation is more complex, as will become clear below.

The methods for neighborhood selection considered here, are based on the column-wise methodology of \cite{meinshausen2006high}. We apply this methodology (under uncertainty) to some known selection methods from the literature, thereby, adjusting these methods for the additional uncertainty. The selection methods considered here are (i) the graphical Lasso of \cite{meinshausen2006high}, (ii) a class of Dantzig-type selectors, that includes the Dantzig selector of \cite{candes2007dantzig}, and (iii) the matrix uncertainty selector of \cite{rosenbaum2010sparse}. This will lead to `graphical' versions of the respective procedures. The graphical Dantzig selector already has been studied in \cite{yuan2010high}, but without the additional uncertainty we are facing here. This leads to novel selection methodologies for which we derive statistical guarantees. We also present numerical studies to illustrate their finite sample performance.

More details on our latent variable block model is discussed in Section~\ref{sec: lbm}. Thereby we also introduce some basic notation. Section~\ref{sec: est} introduces our neighborhood selection method-ologies, and presents results on their large sample performance. Tuning parameter selection is also discussed there. Numerical studies are presented in Section~\ref{sec: num}, and the proofs of our main results are in Appendix~\ref{sec: pfs}..

\section{Some important preliminary facts} \label{sec: lbm}
Let $\bm{\eta}=(\eta_{1},\cdots,\eta_{K})^T$ with $\eta_{k}=\log(p_{kk}/(1-p_{kk}))$ be the vector of log odds of the within-block  connection probabilities, and let $\bm{X}_{K\times L}=(\bm{x}_{1},\cdots,\bm{x}_{K})^T$ be the design matrix. Our latent variable block model \eqref{model: level 2} - \eqref{model: level1_error} says that $\bm{\eta}\sim N(\bm{X}\bm{\beta},\bm{\Sigma})$. The dependence among the $\eta_k,$ encoded in $\bm{\Sigma},$ is propagated to the $p_{kk}$. Let $\bm{\Sigma}^{-1} = \bm{D} = (d_{kl})_{1 \le k,l \le K}$, then the following fact holds.
\begin{fact}
	Under \eqref{model: level 2} - \eqref{model: level1_error}, we have
	$d_{kl}=0,$ if and only if, $p_{k,k}$ is independent of $p_{l,l}$ given the other variables $\bm{p}_{-(k,l)}=\{p_{i,j}: (i, j)\in\wb{K}\times\wb{K}\backslash\{(k, k), (l, l)\}, i \leq j\}$, or just given $\{p_{i,i}: i\in\wb{K}\backslash\{k, l\}\}$.
\end{fact}
\noindent
In other words, if
\begin{align*}
	E=\{(k,l): d_{kl}\neq 0, k\neq l\}
\end{align*}
denotes the edge set of the graph corresponding to $\bm{\eta}$, then, under our latent variable block model, $(k,l)\notin E$ if and only if $p_{k,k}$  is conditionally independent with $p_{l,l}$ given the other variables $\{p_{ii}:1\leq i \leq K, i\notin\{k,l\}\}$. Identifying nonzero elements in $\bm{D}$ thus will reveal the conditional dependence structure of the blocks in our underlying network.

We will use the relative number of edges within each block, as estimates for the unobserved values $p_{kk},k=1,\ldots,K$. Let $S_k=\sum_{\bm{z}[i]=\bm{z}[j]=k}Y_{ij},\,k=1,\cdots, K,$ denote the total number of edges in the $K$ blocks.
\begin{fact} Under \eqref{model: level 2} - \eqref{model: level1_error}, we have
	%
	$${\rm sign}(\sigma_{kl})={\rm sign}\left({\rm Cov}(S_k,S_l)\right).$$
	%
\end{fact}
\noindent
For proofs of the two facts see \cite{paulo2012asymptotics} (page 13, Theorem 1.35) and \cite{liu2009nonparanormal} (Section 3, Lemma 2), respectively.

\section{Neighborhood selection} \label{sec: est}
Here we discuss the identification of the nonzero elements in $\bm{D}$. We first assume that \eqref{model: level 2} - \eqref{model: level1_error} holds with a known $\bm{\beta}$, and we write $\bm{\mu} = (\mu_1,\cdots,\mu_K)^T =\bm{X}\bm{\beta}$. We also assume that $0 < p_{i,j} < 1$ for all $i, j \in\wb{K}$. Let $\bm{Y}^{(t)}, t = 1,\ldots,n$ denote $n$ iid observed networks with corresponding independent unobserved random vectors $\bm{p}^{(t)}, t = 1,\ldots,n$ following our model. Let $\mathcal{A}_1,\cdots,\mathcal{A}_K$ denote the $K$ blocks of the networks $\bm{Y}^{(t)}$ and $\mathcal{V}=\{1,\cdots,N\}$ be the node set.
Assume $\mathcal{A}_k$ and $\mathcal{A}_l$ are mutually exclusive for $k\neq l$ so that $\bigcup_{k=1}^K\mathcal{A}_k=\mathcal{V}$. The number of possible edges within each block is $m_k=|\mathcal{A}_k|(|\mathcal{A}_k|-1)/2$ for $k=1,\cdots,K$, and the number of possible edges between block $k$ and block $l$ is then $|\mathcal{A}_k||\mathcal{A}_l|$ for $1\leq k\neq l\leq K$. We would like to point out again that the block membership variable $\bm{z}$ is assumed to be known.

\subsection{Controlling the estimation error}

Given a network $\bm{Y}^{(t)}$, let $S^{(t)}_k=\sum_{\bm{z}[i]=\bm{z}[j]=k}Y_{ij}^{(t)}, k = 1,\ldots,K, t = 1,\ldots,n$ denote the number of edges within block $k$ in network $t$. Natural estimates of $p_{kk}^{(t)}$ and $\eta_k^{(t)}$ are
\begin{align} \label{def: hatPeta}
	\widetilde{p}_{k,k}^{(t)}=\frac{S_k^{(t)}}{m_k}
	\quad\text{ and}\quad
	\widetilde{\eta}_k^{(t)}=\log\Bigg(\frac{\widetilde{p}_{k,k}^{(t)}}{1-\widetilde{p}_{k,k}^{(t)}}\Bigg)
\end{align}
respectively.

Let $\widetilde{\bm{\eta}}^{(t)}=(\widetilde{\eta}_1^{(t)},\cdots,\widetilde{\eta}_K^{(t)})^T$,
and let $m_{\min}=\min_{1\leq k\leq K}m_k$ be the minimum number of possible edges within a block, which of course measures the minimum blocksize.
\begin{fact}
	Assume that $K$ is fixed. Then, under \eqref{model: level 2} - \eqref{model: level1_error}, we have for each $t = 1,\cdots, n$,
	%
	\begin{align*}
	\widetilde{\bm{\eta}}^{(t)}\rightarrow N(\bm{X}\bm{\beta}, \bm{\Sigma})\quad \text{in distribution as }\,m_{\min}\rightarrow\infty.
	\end{align*}
	%
\end{fact}
This result tells us that, if we base our edge selection on $\widetilde{\bm{\eta}}^{(t)}$, then, for $m_{\min}$ large, we are {\em close} to a Gaussian model, and thus we can hope that our analysis is similar to that of a Gaussian graphical model. However, the approximation error has to be examined carefully. In order to do that, we first truncate the $\widetilde{p}_{kk}^{(t)}$'s, or, equivalently, the $\widetilde{\eta}_k^{(t)}$. For $T > 0,$ let
\begin{equation*}
	\widehat{\eta}_k^{(t)} =
	\begin{cases}
		-T & \mbox{if  $\widetilde{\eta}_k^{(t)} < -T$}\\
		\widetilde{\eta}_k^{(t)} &\mbox{if  $|\widetilde{\eta}_k^{(t)}| \leq T$} \\
		T & \mbox{if  $\widetilde{\eta}_k^{(t)} > T$}.
	\end{cases}
\end{equation*}
This truncation corresponds to
\begin{equation*}
	\widehat{p}_{k,k}^{(t)} =
	\begin{cases}
		(1+e^T)^{-1} & \mbox{if  $\widetilde{p}_{k,k}^{(t)} < (1+e^T)^{-1}$ }\\
		\widetilde{p}_k^{(t)} &\mbox{if  $(1+e^T)^{-1} \leq \widetilde{p}_{k,k}^{(t)} \leq (1+e^{-T})^{-1}$ } \\
		(1+e^{-T})^{-1} & \mbox{if  $\widetilde{p}_{k,k}^{(t)} > (1+e^{-T})^{-1}$}.
	\end{cases}
\end{equation*}
In what follows, we work with these truncated versions. Note that the dependence on $T$ is not indicated explicitly in this notation.

The magnitude of $m_{\min}$  is important, as it reflects the accuracy of our estimates. This estimation error will crucially enter the performance of the graphical model based inverse covariance estimator. Under the latent variable block model, we have the following concentration result:
\begin{lemma} \label{lemma: net_etaDiff}
	Let $\sigma^2 =\max_{k\in\wb{K}}\sigma_{kk}$ and $\mu_B =\max_{k\in\wb{K}}|\mu_k|$. Then, under \eqref{model: level 2} - \eqref{model: level1_error}, we have, for $\min(L,T) > \mu_B$, and $m_{\min}\ge 16M^2\log (nK)e^{2L}$, that
	%
	\begin{multline} \label{eq: net_etaDiff}
		P\Big(\max_{1 \leq k \leq K \atop 1 \leq t \leq n} |\widehat \eta_k^{(t)} - \eta_k^{(t)}| <8Me^{L}\sqrt{\frac{\log (nK) }{m_{\min}}}\Big)  \\
		\geq
		1 -\sqrt{\frac{2}{\pi}}\frac{nK\sigma}{\min\{L, T\}-\mu_B}\exp\Big(-\frac{(\min\{L,T\}-\mu_B)^2}{2\sigma^2}\Big) - \Big(\frac{1}{nK}\Big)^{2M^2 - 1}.
	\end{multline}
	%
\end{lemma}
\begin{remark}
	Note that the larger $\mu_B$, the larger we need to choose both $T$ and $L$. A large $T$ will cause problems, because the $\widehat{p}_k^{(t)}$ then might be too close to zero or one, causing challenges by definition of $\widehat \eta_k^{(t)}$. A large $L$ makes our approximation less tight. Therefore we will have to control the size of $\mu_B$ (even if $\mu_B$ is known); see assumption A1.6 and B1.5.

	To better understand the bound in \eqref{eq: net_etaDiff}, suppose that the number of blocks, $K$, grows with $n$ such that $K(n)=O(n^\gamma),$ for some $\gamma>0.$ While $K$ is allowed to grow with $n$, we assume that $\sigma^2$ is bounded. If we further choose $0 < \min\{L, T\}-\mu_B = \gamma\log n,$ for some $\gamma>0$, then, there exists $c>0$, such that as $n\rightarrow\infty$,
	%
	\begin{align*}
		\sqrt{\frac{2}{\pi}}\frac{nK\sigma}{\min\{L, T\}-\mu_B}\exp\left(-\frac{(\min\{L,T\}-\mu_B)^2}{2\sigma^2}\right) = O\left(\exp\left(-c(\log n)^2\right)\right).
	\end{align*}
	%
	The last term on the right-hand side of \eqref{eq: net_etaDiff} can be controlled similarly, by choosing $M = \sqrt{\left(1+(c\log n)/(\gamma+1)\right)/2}$.
	With these choices, we obtain an approximation error of  $\max_{1 \leq k \leq K, 1 \leq t \leq n} |\widehat \eta_k^{(t)} - \eta_k^{(t)}| = O(n^{-p})$ by choosing the minimum blocksize large enough
	%
	\begin{align*}
		m_{\min}^{-1} = O\left(n^{-2p}(\log n)^{-2}e^{-2L}\right).
	\end{align*}
	%
\end{remark}

\subsection{Edge selection under uncertainty}  \label{subsec: uncertainty}
In order to identify the nonzero elements in $\bm{D}$, we consider the graphical model in terms of the distribution of $\bm{\eta}$.  Recall that $\bm{\eta} \in {\mathbb R}^K$, where each component of $\bm{\eta}$ belongs to one of the $K$ blocks, thus $\wb{K}=\{1,\cdots,K\}$ are not only the block labels, but also the node set in the underlying graph  corresponding to the joint distribution of the $\bm{\eta}$. Using Gaussianity of $\bm{\eta}$, the set ${\rm ne}_a=\{b\in\wb{K}: d_{ab}\neq 0\}$  is the neighborhood of node $a\in\wb{K}$ of the associated graph. We follow the idea of \cite{meinshausen2006high} to convert the problem into a series of linear regression problems: for each $a\in\wb{K}$,
\begin{align*}
	\eta_a - \mu_a = \sum_{b\in\wb{K}\backslash\{a\}}\theta_b^a(\eta_b - \mu_b) + v_a
\end{align*}
with the residual $v_a$ independent of $\{\eta_b: 1\leq b\neq a\leq K\}$. Let $\bm{\theta}^a = (\theta^a_1, \cdots, \theta^a_K)^T \in\mathbb{R}^K$ with $\theta^a_a = 0$, then the neighborhood can also be written as ${\rm ne}_a=\{b\in\wb{K}: \theta_b^a\neq 0\}$.

\cite{meinshausen2006high} consider the case of $n$ i.i.d. observations of $\bm{\eta}$.
However, under the assumption of our model, we only have observations of $\widehat{\bm{\eta}}=(\widehat{\eta}_1,\cdots, \widehat{\eta}_K)^T$. Under our assumptions, we have available $n$ independent realizations $\widehat{\bm{\eta}}^{(1)}, \cdots, \widehat{\bm{\eta}}^{(n)}$.
Let $\widehat{\bm{H}}=(\widehat{\bm{\eta}}^{(1)},\cdots,\widehat{\bm{\eta}}^{(n)})^T$ be the $n\times K$-matrix with columns $\widehat{\bm{\eta}}_a = (\widehat \eta^{(1)}_a,\ldots, \widehat\eta_a^{(n)})^T$, $a \in \wb{K}$.
Similarly denote by $\bm{H} = (\bm{\eta}^{(1)},\ldots,\bm{\eta}^{(n)})^T$ the $n\times K$-matrix whose rows are $n$ independent copies of $\bm{\eta}$. Its column $\bm{\eta}_a, a \in \wb{K}$ are vectors of $n$ independent observations of $\eta_a$. That is, we can also write $\widehat{\bm{H}} = (\widehat{\bm{\eta}}_1,\cdots,\widehat{\bm{\eta}}_K)$ and $\bm{H}=(\bm{\eta}_1,\cdots,\bm{\eta}_K)$.  With this notation, for all $a\in\wb{K}$,
\begin{align} \label{eq: etaV}
	\bm{\eta_a} - \mu_a\bm{1}_n=\sum_{b\in\wb{K}}\theta_b^a(\bm{\eta_b} - \mu_b\bm{1}_n)+\bm{v}_a.
\end{align}
Let $\bm{R}=\widehat{\bm{H}}-\bm{H}$. The new matrix model can be written as
\begin{align} \label{model: matrix}
	&(\widehat{\bm{H}} - \bm{1}_n{\bm{\mu}}^T) = (\bm{H} - \bm{1}_n{\bm{\mu}}^T)+\bm{R} \\
	&(\bm{\eta}^{(t)} - \bm{\mu})\sim N(\bm{0},\bm{\Sigma}) \quad\text{i.i.d. for }t=1,\cdots,n.
\end{align}
Moreover, for each $a\in\wb{K}$, let $\bm{H}_{-a}=\{\bm{\eta}_b: b\in\wb{K}\backslash\{a\}\}$, $\widehat{\bm{H}}_{-a} = \{\widehat{\bm{\eta}}_b: b\in\wb{K}\backslash\{a\} \}$ and  $\bm{\theta}^a_{-a}=(\theta^a_1,\cdots,\theta^a_{a-1},\theta^a_{a+1},\cdots,\theta^a_K)^T$, $\bm{\mu}_{-a} = (\mu_1,\cdots,\mu_{a-1}, \mu_{a+1}, \cdots,\mu_K)^T$. We can write the above model as
\begin{equation} \label{model: mus}
	\begin{aligned}
		\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n = (\bm{H}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T)\bm{\theta}^a_{-a}+\bm{\xi}_a\hspace{7.5mm}\qquad \\
		(\widehat{\bm{H}}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T) = (\bm{H}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T)+\bm{R}_{-a},\qquad
	\end{aligned}
\end{equation}
where $\bm{\xi}_a=\bm{v}_a+(\widehat{\bm{\eta}}_a-\bm{\eta}_a)$ and $\bm{R}_{-a}=\widehat{\bm{H}}_{-a}-\bm{H}_{-a}$.
Note that \eqref{model: mus} has a similar structure as the model considered by \cite{rosenbaum2010sparse}. The important difference is that in our situation, we do not have independence of $\bm{\xi}_a$ and $\bm{R}_{-a}$.

\subsection{Edge selection under uncertainty using the Lasso}
\label{subsec: lasso}
As in \cite{meinshausen2006high}, we define our Lasso estimate $\widehat{\bm{\theta}}^{a,\lambda,lasso}$ of $\bm{\theta}^a$ as
\begin{align} \label{def: thetaLasso}
	\widehat{\bm{\theta}}^{a,\lambda,lasso}=\mathop{\arg\min}\limits_{\bm{\theta}: \theta_a=0}\big(n^{-1}\|(\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n) - (\widehat{\bm{H}} - \bm{1}_n\bm{\mu}^T)\bm{\theta}\|_2^2+\lambda\|\bm{\theta}\|_1\big).
\end{align}
The corresponding neighborhood estimate is
\begin{align*}
	\widehat{{\rm ne}}_a^{\lambda, lasso} = \big\{b\in\wb{K}: \widehat{\theta}^{a, \lambda, lasso}_b\neq0\big\};
\end{align*}
and the full edge set can be estimated by
\begin{align*}
	\widehat{E}^{\lambda,\wedge,lasso} = \big\{(a,b): a\in\widehat{{\rm ne}}_b^{\lambda,lasso} \text{ and } b\in\widehat{{\rm ne}}_a^{\lambda,lasso}\big\}
\end{align*}
or
\begin{align*}
	\widehat{E}^{\lambda,\vee, lasso} = \big\{(a,b): a\in\widehat{{\rm ne}}_b^{\lambda,lasso}\text{ or }b\in\widehat{{\rm ne}}_a^{\lambda,lasso}\big\}.
\end{align*}

In order to formulate statistical guarantees for the behavior of these estimates, we need the following assumptions. On top of the assumptions from \cite{meinshausen2006high}, which are assumptions A1.1 - A1.5, we need further assumption on the underlying network.
\begin{itemize}
	\item[A1] \emph{Assumptions on the underlying Gaussian graph}
	\begin{enumerate}
		\item \emph{High-dimensionality}: There exists some $\gamma>0$ so that
		$K(n)=O(n^\gamma)$ for $n\rightarrow\infty$.

		\item \emph{Nonsingularity}: For all $ a\in\wb{K}$ and $n\in\mathbb{N}$, ${\rm Var}(\eta_a)=1$ and there exists $\upsilon^2>0$ so that
		\begin{align*}
			{\rm Var}(\eta_a|\bm{\eta}_{\wb{K}\backslash\{a\}})
			\geq \upsilon^2.
		\end{align*}

		\item \emph{Sparsity}
		\begin{enumerate}
			\item There exists some $0\leq\kappa<1$ so that $
			\max_{a\in\wb{K}}\left|{{\rm ne}}_a\right|=O(n^\kappa)$ for $n\rightarrow\infty$.
			\item There exists some $\vartheta<\infty$ so that for all neighboring nodes $a,b\in\wb{K}$ and all $n\in\mathbb{N}$,
			\begin{align*}
				\|\bm{\theta}^{a,{{\rm ne}}_b\backslash\{a\}}\|_1\leq\vartheta.
			\end{align*}
		\end{enumerate}

		\item \emph{Magnitude of partial correlations}: There exist a constant $c>0$ and some $1\geq\xi>\kappa$, so that for all $(a,b)\in E$,
		\begin{align*}
			|\pi_{ab}|\geq cn^{-(1-\xi)/2},
		\end{align*}
		where $\pi_{ab}$ is the partial correlation between $\eta_a$ and $\eta_b$.

		\item \emph{Neighborhood stability}: There exists some $\varrho<1$ so that for all  $ a,b\in\wb{K}$ with $b\notin{{\rm ne}}_a$,
		\begin{align*}
			|S_a(b)|<\varrho
		\end{align*}
		where
		\begin{align*}
			S_a(b)=\sum_{k\in{{\rm ne}}_a}{\rm sign}\left(\theta_k^{a,{{\rm ne}}_a}\right)\theta_k^{b,{{\rm ne}}_a}.
		\end{align*}

		\item \emph{Asymptotic upper bound on the mean:} $\mu_B(n)=o(\log n)$ for $n\rightarrow\infty$.
	\end{enumerate}

	\item[A2] \emph{Block size of networks}:
	There exists constants $c>0$ and $n_0,$ such that
	\begin{align*}
		m_{\min}(n) \ge c\cdot n^\nu
		\quad
		\text{for }n\geq n_0,
	\end{align*}
	where $\nu>\max\{4-4\xi,2-2\xi+2\kappa\}.$
\end{itemize}
The following theorem shows that, for proper choice of $\lambda=\lambda_n$, our selection procedure finds the correct neighborhoods with high probability, provided $n$ is large enough.
\begin{theorem} \label{corollary: net_mb0_asy}
	Let assumptions A1 and  A2 hold, and assume $\bm{\beta}$ to be known. Let $\epsilon$ be such that
	\begin{align*}
		0 < \max\big\{\kappa, \tfrac{4-\xi-\nu}{3},\tfrac{2+2\kappa-\nu}{2} \big\} < \epsilon <\xi.
	\end{align*}
	If, for some $d_T, d_\lambda >0,$ we have $T_n \sim d_T\log n$ and $\lambda_n\sim d_\lambda n^{-(1-\epsilon)/2}$ \footnote{For two sequence $\{a_n\}$, $\{b_n\}$ of real numbers, we write $a_n \sim b_n$ for $\frac{a_n}{b_n}\rightarrow c$ for some $0 < c <\infty$. },
	respectively, then there exists a constant $c>0,$ such that
	\begin{align*}
		P(\widehat{E}^{\lambda,lasso} = E)
		= 1 - O\left(\exp\left(-c(\log n)^2\right)\right)
		\quad
		\text{as }n\rightarrow\infty.
	\end{align*}
\end{theorem}
\begin{remark}
	Assumption A2 says that the rate of increase of the minimum block size, which behaves like $\sqrt{m_{\min}},$ depends on the neighborhood size in  our graphical model, and on the magnitude of the partial correlations in the graphical model. Roughly speaking, large neighborhoods (large $\kappa$), and small partial correlations (small $\xi$), both require a large minimum block size (large $\nu$), which appears reasonable. The choice of a proper penalty parameter $\lambda_n$ also depends on these two parameters.
\end{remark}

\subsection{Edge selection with a class of Dantzig-type selectors under uncertainty} \label{subsec: ds}
In this section, we propose a novel class of Dantzig-type selectors that are iterated over all $a\in\wb{K}$. For a linear model as in \eqref{eq: etaV}, i.e. for fixed $a$, \cite{candes2007dantzig} introduced the Dantzig selector as a solution to the convex problem
\begin{multline*}
	\quad
	\min\Big\{\|\bm{\theta}\|_1:\bm{\theta}\in\mathbb{R}^{K(n)}, \theta_a=0\text{ and } \\
	\Big|\frac{1}{n}\big(\bm{H}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T\big)^T\big((\bm{\eta}_a - \mu_a\bm{1}_n) - (\bm{H} - \bm{1}_n{\bm{\mu}}^T)\bm{\theta}\big)\Big|_\infty \leq \lambda\Big\},
	\quad
\end{multline*}
where $\lambda\geq0$ is a  tuning parameter, and for a matrix $A =(a_{ij})$, $|\cdot|_\infty = \max_{ij}|a_{ij}|.$  Under our model, we define the Dantzig selector as a solution of the minimization problem
\begin{multline} \label{def: thetaDs}
	\quad
	\min\Big\{\|\bm{\theta}\|_1:\bm{\theta}\in\mathbb{R}^{K(n)}, \theta_a=0\text{ and } \\
	\Big|\frac{1}{n}\big(\widehat{\bm{H}}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T\big)^T\big((\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n) - (\widehat{\bm{H}} - \bm{1}_n\bm{\mu}^T)\bm{\theta}\big)\Big|_\infty \leq \lambda\Big\}
	\quad
\end{multline}
with $\lambda\geq0$. Moreover, when considering \eqref{model: mus}, the idea of matrix uncertainty selector (MU-selector) comes into our mind. In our setting, we define an MU-selector, a generalization of the Dantzig selector under matrix uncertainty, as a solution of the minimization problem
\begin{multline} \label{def: thetaMus}
	\quad \min\Big\{\|\bm{\theta}\|_1:\bm{\theta}\in\mathbb{R}^{K(n)}, \theta_a=0\text{ and} \\
	\Big|\frac{1}{n}\big(\widehat{\bm{H}}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T\big)^T\big((\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n)-(\widehat{\bm{H}} - \bm{1}_n\bm{\mu}^T)\bm{\theta}\big)\Big|_\infty \leq \mu\|\bm{\theta}\|_1+\lambda\Big\}
\end{multline}
with tuning parameters $\mu\geq0$ and $\lambda\geq0$. Note that our MU-selector deals with matrix uncertainty directly, rather than replacing $\bm{H}$ by $\widehat{\bm{H}}$ in the optimization equations like the Lasso or the Dantzig selector. What we mean by this is that our MU-selector is based on the structural equation \eqref{model: mus}, while both Lasso-based estimator and Dantzig selector are based on the linear model \eqref{eq: etaV} with the unknown $\bm{\eta}$'s simply replaced by their estimators.

Now we consider a class of Dantzig-type selectors, which can be considered as generalizations of the Dantzig selector and the MU-selector.  For each $a\in\wb{K}$, let the Dantzig-type selector $\widetilde{\bm{\theta}}^{a,\lambda,ds}$ be a solution of the optimization problem
\begin{multline} \label{def: thetaG}
	\quad \min\Big\{\|\bm{\theta}\|_1: \bm{\theta}\in\mathbb{R}^{K(n)}: \theta_a=0\text{ and } \\
	\Big|\frac{1}{n}(\widehat{\bm{H}}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T)^T\big((\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n)-(\widehat{\bm{H}} - \bm{1}_n\bm{\mu}^T)\bm{\theta}\big)\Big|_\infty
	\leq \lambda_{a,n}(\|\bm{\theta}\|_1)\Big\},
\end{multline}
where for each $n\in\mathbb{N}$, $\{\lambda_{a,n}(\cdot): a\in\wb{K}\}$ is a set of functions such that
\begin{itemize}
	\item For each $n\in\mathbb{N}$ and $a\in\wb{K}$, $\lambda_{a,n}(\cdot)$ is an increasing function.

	\item For all $n\in\mathbb{N}$, $\min_{a\in\wb{K}}\lambda_{a,n}(\cdot)$ is lower bounded by some constant $\lambda_n \geq 0$, i.e, for all $n\in\mathbb{N}$, there exists some $\lambda_n >0$ so that
	\begin{align*}
		\min_{ a\in\wb{K}}\min_{\bm{\theta}\in\mathbb{R}^{K}: \theta_a = 0}\lambda_{a,n}(\|\bm{\theta}\|_1) \geq \lambda_n.
	\end{align*}

	\item $\max_{a\in\wb{K}}\lambda_{a,n}(\|\bm{\theta}^a\|_1) = o(n^{-\frac{1-\xi}{2}}(\log n)^{-1})$, i.e, there exist $u_n = o(1)$ and $n_0\in\mathbb{N}$, so that, for all $n\geq n_0$,
	%
	\begin{align*}
		\lambda_{a,n}(\|\bm{\theta}^a\|_1) \leq u_nn^{-\frac{1-\xi}{2}}(\log n)^{-1}, \text{ for all }a\in\wb{K}.
	\end{align*}
	%
\end{itemize}
The Dantzig-type selector $\widetilde{\bm{\theta}}^{a,\lambda,ds}$ always exists, because the LSE $\widehat{\bm{\theta}}^a$ defined as $\widehat{\bm{\theta}}^a_{-a} = (\widehat{\bm{H}}_{-a} -\bm{1}_n\bm{\mu}_{-a}^T )^{+}(\hat{\bm{\eta}}_a - \mu_a\bm{1}_n)$ and $\widehat{\theta}^a_{a} = 0$ belongs to the feasible set $\Theta_a$, where
\begin{multline*}
	\Theta_a=\Big\{\bm{\theta}\in\mathbb{R}^{K(n)}: \theta_a=0\text{ and } \\
	\Big|\frac{1}{n}\big(\widehat{\bm{H}}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T\big)^T\big((\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n)-(\widehat{\bm{H}} - \bm{1}_n\bm{\mu}^T)\bm{\theta}\big)\Big|_\infty \leq \lambda_{a,n}(\|\bm{\theta}\|_1)\Big\}
\end{multline*}
for any $\lambda_{a,n}(\|\bm{\theta}\|_1)\geq0$. It may not be unique, however.  We will show that, similar to \cite{candes2007dantzig} and \cite{rosenbaum2010sparse}, under certain conditions, for large $n$, there exists a constant $t > 0$ such that the $l_\infty$-norm of the difference between the Dantzig-type selector $\widetilde{\bm{\theta}}^{a,\lambda,ds}$ and the population quantity $\bm{\theta}^a,$ can be bounded by $t\lambda_{a,n}(\|\widetilde{\bm{\theta}}^{a,\lambda,ds}\|_1)$ for all $a\in\wb{K}$ with large probability, where $t_n$ can be a constant large enough or of order $\log n$.  However, in general, sparseness cannot be guaranteed. This already has been observed in \cite{rosenbaum2010sparse}. Therefore, we consider a thresholded version of the Dantzig-type selector, which can also significantly improve the accuracy of the estimation of the sign.
Let $\widehat{\bm{\theta}}^{a,\lambda,ds}\in\mathbb{R}^{K(n)}$ be defined as
\begin{align} \label{def: thetaPd}
	\widehat{\theta}^{a,\lambda,ds}_b=
	\begin{cases}
		0, & \quad b=a \\
		\widetilde{\theta}^{a,\lambda,ds}_bI\big(|\widetilde{\theta}^{a,\lambda,ds}_b| > t_n\lambda_{a,n}(\|\widetilde{\bm{\theta}}^{a,\lambda,ds}\|_1)\big), & \quad
		b\in\wb{K}\backslash\{a\}, \\
	\end{cases}
\end{align}
where $I(\cdot)$ is the indicator function, and $t_n$ is a sequence that satisfies $t_n^{-1} = o(1)$ and $t_n = O(\log n)$.  The corresponding neighborhood selector is, for all $a\in\wb{K},$  defined as $\widehat{{\rm ne}}_a^{\lambda,ds}=\{b\in\wb{K}: \widehat{\theta}_b^{a,\lambda,ds}\neq 0\},$
and the corresponding full edge selector is
\begin{align*}
	\widehat{E}^{\lambda,\wedge,ds}=\big\{(a,b): a\in\widehat{{\rm ne}}_b^{\lambda,ds} \text{ and } b\in\widehat{{\rm ne}}_a^{\lambda,ds}\big\}
\end{align*}
or
\begin{align*}
	\widehat{E}^{\lambda,\vee,ds}=\big\{(a,b): a\in\widehat{{\rm ne}}_b^{\lambda,ds} \text{ or } b\in\widehat{{\rm ne}}_a^{\lambda,ds}\big\}.
\end{align*}

Similar to the Section~\ref{subsec: lasso}, in order to derive some consistency properties, we need assumption about the underlying Gaussian graph (B1), and the minimum block size in the underlying network (B2).
\begin{itemize}
	\item[B1] \emph{Assumptions on the underlying Gaussian graph}
	\begin{enumerate}
		\item \emph{Dimensionality}: There exists $\gamma>0$ such that
		$K(n)=O(n^\gamma)$ as $n\rightarrow\infty$.

		\item \emph{Nonsingularity}: For all $ a\in\wb{K}$ and $n\in\mathbb{N}$, ${\rm Var}(\eta_a)=1$ and there exists $\upsilon^2>0$ so that
		\begin{align*}
		{\rm Var}(\eta_a|\bm{\eta}_{\wb{K}\backslash\{a\}})\geq \upsilon^2.
		\end{align*}

		\item \emph{Sparsity}
		\begin{enumerate}
			\item There exists
			$0\leq\kappa<1/2$, so that $
			\max_{a\in\wb{K}}\left|{{\rm ne}}_a\right|=O(n^\kappa)$, as $n\rightarrow\infty$.
			\item $\|\bm{\Sigma}^{-1}\|_\infty = O(1)$ as $n\rightarrow\infty$.
		\end{enumerate}

		\item \emph{Magnitude of partial correlations}: There exist a constant $c>0$ and $1\geq\xi>\kappa$, so that, for all $(a,b)\in E$, $|\pi_{ab}|\geq cn^{-(1-\xi)/2}$.

		\item \emph{Asymptotic upper bound on the mean:} $\mu_B(n)=o(\log n)$ for $n\rightarrow\infty$.
	\end{enumerate}

	\item[B2] \emph{Block size of networks}: $m_{\min}^{-1}(n)=O(n^{-\nu})$ with some $\nu>1-\xi+2\kappa$ for $n\rightarrow\infty$.
\end{itemize}
Here, the assumption on $m_{\min}$ (assumption B2) is weaker than that assumed for the Lasso-based estimator (assumption A2). Similar remarks as given for A2 also apply to B2 (see Remark right below Theorem~\ref{corollary: net_mb0_asy}).

Assumptions A1 and B1 are similar but not equivalent: A1.1 and B1.1, A1.2 and B1.2, A1.4 and B1.4 respectively, are exactly the same. B1.2(a)  is stronger than A1.3.(a), indicating the underlying graph should be even sparser than the graph in Section~\ref{subsec: lasso}; assumption B1 does not have an analog to A1.3.(b) and A1.5.
\begin{theorem} \label{corollary: net_ds0_asy}
	Let assumptions B1 and B2 hold, and assume $\bm{\beta}$ is known. Let $\epsilon>0$ be such that $\xi>\epsilon>1+2\kappa-\nu$. If $T_n \sim d_T\log n$ with some $d_T>0$, and  $\lambda_n^{-1} = O(n^{\frac{1-\epsilon}{2}})$, there exists $c>0$, so that
	\begin{align*}
		P(\widehat{E}^{\lambda,ds} = E)
		= 1 - O\left(\exp\left(-c(\log n)^2\right)\right)
		\quad
		\text{as }n\rightarrow\infty.
	\end{align*}
\end{theorem}
\begin{remark}
	The choice of proper $\lambda_{a,n}(\cdot)$ depends on the three parameters $\xi, \kappa$ and $\nu$. However, even the best scenario does not allow for the order $\lambda_p \sim \sqrt{\frac{\log p}{n}},$ which often can be found in the literature. This stems from the fact that we have to deal with an additional estimation error (coming in through the estimation of $\boldsymbol \eta$).
\end{remark}

\subsection{Extension} \label{subsec: extension}
Here we consider the case of an unknown coefficient vector $\bm{\beta},$ or unknown mean $\bm{\mu} = \bm{X}\bm{\beta}$.
Recall that $\bm{\eta}^{(t)}\sim N(\bm{\mu},\bm{\Sigma}),$ $t=1,\cdots,n$ are i.i.d. Given $\{\bm{\eta}^{(t)}: t=1,\cdots,n\}$, a natural way to estimate $\bm{\mu}$ is via the MLE $\widebar{\bm{\eta}}=\frac{1}{n}\sum_{t=1}^{n}\bm{\eta}^{(t)}$. Recall, however, that we only have estimates $\widehat{\bm{\eta}}^{(t)}, t=1,\cdots,n$, available. Using the estimates $\widehat{\bm{\eta}}^{(t)}$, we estimate the underlying mean $\bm{\mu}$ by  $\widebar{\widehat{\bm{\eta}}}=\frac{1}{n}\sum_{t=1}^{n}\widehat{\bm{\eta}}^{(t)}$. Moreover, we can estimate $\bm{\beta}$ via $\widehat{\bm{\beta}}=\bm{X}^{+}\widebar{\widehat{\bm{\eta}}}$, where $\bm{X}^+$ is the Moore-Penrose pseudoinverse of $\bm{X}$ (when ${\rm rank}(\bm{X})=L$, $\bm{X}^+=(\bm{X}^T\bm{X})^{-1}\bm{X}^T$). In order to derive consistency properties for $\widehat{\bm{\beta}}$, assumptions on the design matrix are needed.
Theorem~\ref{thm: net_beta} below states asymptotic properties of the estimators.
\begin{theorem} \label{thm: net_beta}
	Let assumptions A1.1 (or B1.1) and A1.6 (or B1.5) hold. If $m^{-1}_{\min}=O(n^{-\nu})$ for some $\nu>0$, then, for any $b<\min\{1,\nu\}$, and fixed $\delta>0$, there exists some $c>0$ so that
	\begin{align*}
		P\big(n^{\tfrac{b-\gamma}{2}}\|\widebar{\widehat{\bm{\eta}}}-\bm{\mu}\|_2>\delta\big)
		=
		O\left(\exp\left(-c(\log n)^2\right)\right)
		\quad
		\text{as }n\rightarrow\infty.
	\end{align*}
	If, moreover, the design matrix is of full rank and the singular value of $\bm{X}$ is asymptotically upper bounded, that is, ${\rm rank}(\bm{X})=L$ and $\sigma_{\max}(\bm{X})=O(1)$, then there exists $c>0,$ so that
	\begin{align*}
		P\big(n^{\tfrac{b-\gamma}{2}}\|\widehat{\bm{\beta}}-\bm{\beta}\|_2>\delta\big)=O\left(\exp\left(-c(\log n)^2\right)\right)
		\quad
		\text{as }n\rightarrow\infty.
	\end{align*}
\end{theorem}
Next we consider the estimation of the edge set $E$ based on $\bm{D}=\bm{\Sigma}^{-1}$. We write $\bm{\eta}-\bm{\mu}\sim N(\bm{0},\bm{\Sigma})$ and consider $(\widehat{\bm{\eta}}^{(t)}-\widebar{\widehat{\bm{\eta}}})_{t=1,\cdots,n}$ as the observations. We estimate the edge set in the same way as described in Section~\ref{subsec: lasso}, but replace $\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n$ by $\widehat{\bm{\eta}}_a-\widebar{\widehat{\eta}}_a\bm{1}_n$ and replace $\widehat{\bm{H}}$ by $\widehat{\bm{H}}-\bm{1}_n\widebar{\widehat{\bm{\eta}}}^{T}$ in \eqref{def: thetaLasso}, where $\widebar{\widehat{\eta}}_a=\frac{1}{n}\sum_{t=1}^n\widehat{\eta}_a^{(t)}$ and $\widebar{\widehat{\bm{\eta}}}$ is as above.
The following consistency result parallels Theorem~\ref{corollary: net_mb0_asy} and Theorem~\ref{corollary: net_ds0_asy}, but stronger assumption are needed to control the additional estimate error.
\begin{corollary} \label{corollary: net_mb_asy}
	Let assumptions A1 - A2 hold with $\xi>3/4$, and let $\epsilon$ be such that
	\begin{align*}
		\max\{\kappa+1/2, \tfrac{3-\xi}{3}, \tfrac{4-\xi-\nu}{3},\tfrac{2+2\kappa-\nu}{2}\} < \epsilon < \xi.
	\end{align*}
	Suppose that $T_n \sim d_T\log n$, for $d_T>0,$ and that the penalty parameter satisfies $\lambda_n\sim d_\lambda n^{-(1-\epsilon)/2}$ for some $d_\lambda>0$. Then, there exists $c>0,$ so that
	\begin{align*}
		P(\widehat{E}^{\lambda,lasso}=E)=1- O\left(\exp\left(-c(\log n)^2\right)\right)
		\quad
		\text{as }n\rightarrow\infty.
	\end{align*}
\end{corollary}
\begin{corollary} \label{corollary: net_ds_asy}
	Let assumptions B1 - B2 hold with $\xi>2\kappa$. Let $\epsilon$ be such that  $\xi>\epsilon>\max\{2\kappa, 2\kappa+1-\nu\}$. If $T_n \sim d_T\log n,$ for some $d_T>0$, and  $\lambda_n^{-1} = O(n^{\frac{1-\epsilon}{2}})$, there exists $c>0$ so that
	\begin{align*}
		P(\widehat{E}^{\lambda,ds}=E)=1- O\left(\exp\left(-c(\log n)^2\right)\right)
		\quad
		\text{as }n\rightarrow\infty.
	\end{align*}
\end{corollary}
\begin{example} \label{eg: all1}
	Let $\kappa=0$ and $\xi=1$, that is, the number of blocks is finite, and the partial correlations are lower bounded for the graphcial model. If, in addition, for some $\nu > 0,$ $m_{\min}^{-1}(n)=O(n^{\nu})$ as $n\rightarrow\infty$, then Corollaries~\ref{corollary: net_mb_asy} and \ref{corollary: net_ds_asy}, respectively, apply in the following scenarios:
	%
	\begin{itemize}
		\item \textbf{The Lasso}: If assumption A1 - A2 hold: Choose the tuning parameter $\lambda_n\sim dn^{-(1-\epsilon)/2}$ with any $\epsilon$ satisfying $1>\epsilon>\max\{0,1-\nu/2\}$ in case $\bm{\mu}$ is known, and $\epsilon$ satisfying $1>\epsilon>\max\{2/3,1-\nu/2\}$ for $\bm{\mu}$ unknown.

		\item \textbf{The Dantzig-type selector}: If assumptions B1 - B2 hold, whether $\bm{\mu}$ is known or unknown,
		choose $\max_{a\in\wb{K}}\lambda_{a,n}^{-1}(\bm{\theta}^a) = O(n^{\frac{1-\epsilon}{2}})$ with any positive $\epsilon$ satisfying $1>\epsilon>1-\nu$. In particular, for \\[-15pt]
		%
		\begin{itemize}[label=\textasteriskcentered]
			\item \textbf{the Dantzig selector}: $\lambda_{a,n}(\|\bm{\theta}\|_1) = dn^{-\frac{1-\epsilon}{2}}$ for any $d>0$. Problem \eqref{def: thetaG} becomes  \eqref{def: thetaDs} with tuning parameter $\lambda = dn^{-\frac{1-\epsilon}{2}}$.

			\item \textbf{the MU-selector}: $\lambda_{a,n}(\|\bm{\theta}\|_1) = dn^{-\frac{1-\epsilon}{2}}\|\bm{\theta}\|_1+dn^{-\frac{1-\epsilon}{2}}$ for any $d>0$. Problem \eqref{def: thetaG} becomes \eqref{def: thetaMus} with tuning parameter $\mu=dn^{-\frac{1-\epsilon}{2}}$ and $\lambda = dn^{-\frac{1-\epsilon}{2}}$.
		\end{itemize}
		%
	\end{itemize}
\end{example}

\subsection{Selection of penalty parameters in finite samples.} \label{subsec: finite}
The results above only show that consistent edge selection is possible with the Lasso and the Dantzig-type selector in a high-dimensional setting. However, we still have not given a concrete way to choose the penalty parameter for a given data set. In this section, we discuss the choice of tuning parameter for finite $n$ for the following estimation methods:
\begin{itemize}
	\item The Lasso
	\item The Dantzig-type selectors:
	\begin{itemize}
		\item the Dantzig selector: $\lambda(\|\bm{\theta}\|_1) = \lambda$
		\item the MU-selector: $\lambda(\|\bm{\theta}\|_1) = \lambda\|\bm{\theta}\|_1 + \lambda$
	\end{itemize}
\end{itemize}

\cite{meinshausen2006high} discussed a data-driven choice of the penalty parameter of the Lasso for Gaussian random vectors. Our data are not Gaussian, however. Moreover, according to our numerical studies,
the choice suggested by Meinshausen and B\"{u}hlmann tends to result in a very sparse graph, which goes along with a very small type I error. Another natural idea is choosing the penalty parameter via cross-validation. However, \cite{meinshausen2006high} already state that the choice $\lambda_{\rm oracle}$ gives an inconsistent estimator, and $\lambda_{\rm cv}$ is an estimate of $\lambda_{\rm oracle}$. So the cross-validation approach is also not recommended. Instead we here consider the following two-stage procedure: for each $a\in\wb{K}$, let
\begin{equation} \label{def: theta2stage}
	\widehat{\theta}_b^{a, \lambda, \tau} =
	\begin{cases}
		\widetilde{\theta}_b^{a, \lambda}I\big(|\widetilde{\theta}_b^{a, \lambda}|/\max_{k\in\wb{K}\backslash\{a\}}|\widetilde{\theta}_b^{a, \lambda}| > \tau \big) & \mbox{if  $\widetilde{\bm{\theta}}^{a, \lambda} \neq \bm{0}$ }\\

		0 & \mbox{if  $\widetilde{\bm{\theta}}^{a, \lambda} = \bm{0}$ },
	\end{cases}
\end{equation}
where $\widetilde{\bm{\theta}}^{a, \lambda}$ is obtained by solving either \eqref{def: thetaLasso}, \eqref{def: thetaDs} or \eqref{def: thetaMus} with $\mu = \lambda$. Such procedures have also been used in \cite{rosenbaum2010sparse} and \cite{zhou2011high}. However, the use of $\max_{k\in\wb{K}\backslash\{a\}}|\widehat{\theta}_b^{a, \lambda}|$ in the truncation is novel. By using $\max_{k\in\wb{K}\backslash\{a\}}|\widehat{\theta}_b^{a, \lambda}|$, we have $\tau \in [0, 1]$, making the tuning parameter $\tau$ more standardized.
Note that when  $\widehat{\theta}_b^{a, \lambda}$ is a Dantzig-type selector, then, under the assumptions in Section~\ref{sec: est}, and for large $n$, \eqref{def: theta2stage} is equivalent to \eqref{def: thetaPd}.

For the choice of $\lambda$ and $\tau$, we follow a similar idea as in \cite{zhou2011high}, but with some modification: for each $a\in\wb{K}$, we select $\lambda_a$ via cross-validation to minimize the squared error prediction loss for $a$-th regression. After all $\lambda_a$,  $a\in\wb{K}$, are chosen, we select $\tau$ via BIC based on a Gaussian assumption:
\begin{align*}
	{BIC}(\bm{D}) = -2l_n(\bm{D}) + \log(n)\dim(\bm{D}),
\end{align*}
where $l_n(\bm{D})$ is the $n$-sample Gaussian log-likelihood and $\dim(\bm{D}) = $ number of free parameters. Note that we do not have a nice form of the likelihood, so we use the Gaussian likelihood instead.

\section{Simulation study} \label{sec: num}
In this section, we mainly study the finite sample behavior of the three estimation methods mentioned in Section~\ref{subsec: finite}, that is,
\begin{itemize}
	\item the Lasso;
	\item the Dantzig selector;
	\item the MU-selector with $\mu = \lambda$.
\end{itemize}

\subsection{Finite-sample performance as a function of the penalty parameter.} \label{subsec: finite_performance}
Here we consider the methods proposed in Section~\ref{subsec: lasso} and \ref{subsec: ds} with an AR(1) type covariance structure $\bm{\Sigma}_{K\times K} = \{\rho^{|i-j|}\}_{i, j \in \wb{K}}$ with $\rho = 0.2, 0.5$ and $0.8$. In this setup, $d_{ij} = 0$ if and only if $|i - j| > 1$. The minimum blocksize in our simulation is set to be $m_{\min} = 100$. We consider the following choices of the sample size and number of blocks:
\begin{itemize}
	\item $n = 20$ with $K = 15$;
	\item $n = 100$ with $K = 15, 30, 50, 80, 100$ and $150$.
\end{itemize}
We only present the results for $n = 20, K = 15$ and $n=100, K=150.$ The rest of the results can be found in the supplementary material.
Figures~\ref{fig: roc1} - \ref{fig: roc7} show ROC-curves; average error rates (total error, type \uppercase\expandafter{\romannumeral1} error and type \uppercase\expandafter{\romannumeral2} error) as functions of the tuning parameter $\lambda$ are shown in figures~\ref{fig: rate1}  and \ref{fig: rate7}.
The shown curves are color-coded: Lasso: red, Dantzig selector: blue and MU-selector: green. $\lambda_{\rm opt}$ is the tuning parameter corresponding to the total (overall) minimum error rate.
\begin{figure} [H]
	\centering
	\includegraphics[width=1\linewidth]{roc1}
	\vspace*{-1cm}
	\caption{\small ROC curves comparing the three proposed methods for $K = 15$ and $n = 20$}
	\label{fig: roc1}
\end{figure}
\begin{figure} [H]
	\centering
	\includegraphics[width=1\linewidth]{roc7}
	\vspace*{-1cm}
	\caption{\small ROC curves comparing the three proposed methods for $K = 150$ and $n = 100$}
	\label{fig: roc7}
\end{figure}
We can see that the value of $\rho$ is important. The performance of all the three methods improves as $\rho$ grows. This can be understood by the fact that it determines the size of the partial correlations (cf. assumption A1.4).
\begin{figure} [H]
	\centering
	\includegraphics[width=1\linewidth]{rate1}
	\caption{\small Average error rates as functions of $\lambda$ for $K = 15$ and $n = 20$}
	\label{fig: rate1}
\end{figure}
Moreover, when $n \geq K$, and for $\lambda = 0$, all these methods result in estimates with all components being non-zero, which result in type \uppercase\expandafter{\romannumeral1} error rate equal to $1$ and type \uppercase\expandafter{\romannumeral2} error equals $0$, that is, $(1, 1)$ in the ROC curves. However, when $n < K$, and $\lambda = 0$, the feasible set $\Theta_a$ is dimension at least $(K - n)$.  The Dantzig-type selectors minimize the $L_1$-norm of these $\bm{\theta}$'s, which produces some zero terms of the solution; thus, the corresponding type \uppercase\expandafter{\romannumeral1}  error rate will be less than $1$ and the type \uppercase\expandafter{\romannumeral2} error rate might be greater than $0$, that is why the ROC curves of the Dantzig selector and the MU-selector cannot reach $(1, 1)$ for the case $n = 100$ with $K = 150$. However, the solution of the Lasso is not unique, the coordinate decent algorithm could return a solution with all its elements non-zero, resulting $(1, 1)$ in the ROC curves.
\begin{figure} [H]
	\centering
	\includegraphics[width=1\linewidth]{rate7}
	\caption{Average error rates as functions of  $\lambda$ for $K = 150$ and $n = 100$}
	\label{fig: rate7}
\end{figure}

\subsection{Finite-sample performance with data-driven penalty selection} \label{subsec: finite_penalty}
In this section, we  study the three methods for finite-sample setup discussed in Section~\ref{subsec: finite}. In our simulation study, we consider three different covariance models with $K = 30, 100, 200$, $m_{\min} = 45$ and $n = 100, 500, 1000$. Below we only present the case $K=100$. See supplemental material for the other cases.
\begin{itemize}
	\item AR(1): $\bm{\Sigma}_{K\times K} = \{\rho^{|i-j|}\}_{1\leq i, j \leq K}$ with $\rho = 0.7$.

	\item AR(4): $d_{ij} = I(|i - j| = 0) + 0.4\dot I(|i - j| = 1) + 0.2\cdot I(|i - j| = 2) + 0.2\cdot I(|i - j| = 3) + 0.1\cdot I(|i - j| = 4)$.

	\item A random precision matrix model (see \cite{rothman2008sparse}): $\bm{D}_K = \bm{B} + \delta\bm{I}$ with each off-diagonal entry in $\bm{B}$ is generated independently and equals $0.5$ with probability $\alpha$ or $0$ with probability $1 - \alpha$. $\bm{B}$ has zeros on the diagonal, and $\delta$ is chosen so that the condition number of $\bm{D}$ is $K$.
\end{itemize}
As mentioned in Section~\ref{subsec: finite}, we choose $\lambda$ via cross-validation, and $\tau$ based on $BIC(\tau)$. As for the choice of $\tau$, we often encountered the problem of a very flat BIC-function close to the level of the minimum (some $BIC(\tau)$ plots are shown in figure~\ref{fig: bic}). To combat this problem, we use the following strategy in our simulations: if more than half of the $\tau\in[0, 1]$ result in the same BIC, then we choice the third quartile of these $\tau$'s, otherwise, we choose the one resulting the minimum BIC.
\begin{figure} [H]
	\begin{subfigure}[b]{0.5\linewidth}
		\centering
		\includegraphics[width=1\linewidth]{bic1}
	\end{subfigure
	\begin{subfigure}[b]{0.5\linewidth}
		\centering
		\includegraphics[width=1\linewidth]{bic2}

	\end{subfigure}
	\begin{subfigure}[b]{0.5\linewidth}
		\centering
		\includegraphics[width=1\linewidth]{bic3}
	\end{subfigure
	\begin{subfigure}[b]{0.5\linewidth}
		\centering
		\includegraphics[width=1\linewidth]{bic4}
	\end{subfigure}
	\caption{Tuning parameter $\tau$ vs $BIC(\tau)$ plots for the four models with $K = 30$: the red points corresponds to the optimal $\tau$}
	\label{fig: bic}
\end{figure}

Simulation results for AR(1) and AR(4) models are shown in tables~\ref{tab: ar1_k100} and \ref{tab: ar4_k100}, respectively. For the random precision matrix model we consider $\alpha = 0.1$ and $\alpha = 0.5$ (as in \cite{zhou2011high}). The simulation results are shown in tables~\ref{tab: rd0.1_k100} and \ref{tab: rd0.5_k100}, respectively. The tables show averages and SEs of classification errors in \% over 100 replicates for the three proposed methods with both $``\vee"$ (left) and $``\wedge"$ (right).

\begin{table} [H]
	\caption{AR(1) model with $K = 100$}
	\begin{subtable}{1\textwidth}
		\sisetup{table-format=-1.2}
		\centering
		\caption{$n = 100$}
		\begin{tabular}{c|c c c }
			\hline
			Ave (SE)
			& Total (\%) &  Type \uppercase\expandafter{\romannumeral1} (\%)  & Type \uppercase\expandafter{\romannumeral2} (\%) \\
			\hline
			{Lasso}&0.90(0.45); 1.16(0.44)&0.91(0.47); 1.18(0.45)&0.39(0.68); 0.06(0.28)\\
			{Dantzig} &0.84(0.42); 1.26(0.44)&0.85(0.43); 1.28(0.45)&0.38(0.68); 0.12(0.36)
			\\
			{MU} &0.93(0.49); 1.22(0.43)&0.94(0.50); 1.24(0.44)&0.23(0.53); 0.15(0.44)\\
			\hline
		\end{tabular}
	\end{subtable}

	\bigskip

	\begin{subtable}{1\textwidth}
		\sisetup{table-format=-1.2}
		\centering
		\caption{$n = 500$}
		\begin{tabular}{c|c c c }
			\hline
			Ave (SE)
			& Total (\%) &  Type \uppercase\expandafter{\romannumeral1} (\%)  & Type \uppercase\expandafter{\romannumeral2} (\%) \\
			\hline
			{Lasso}&0.566(0.220); 0.716(0.231)&0.577(0.224); 0.730(0.236)&0(0)\\
			{Dantzig} &0.601(0.208); 0.700(0.255)&0.613(0.212); 0.714(0.260)&0(0)\\
			{MU} &0.574(0.203); 0.676(0.238)&0.586(0.207); 0.690(0.243)&0(0)\\
			\hline
		\end{tabular}
	\end{subtable}

	\bigskip

	\begin{subtable}{1\textwidth}
		\sisetup{table-format=-1.2}
		\centering
		\caption{$n = 1000$}
		\begin{tabular}{c|c c c }
			\hline
			Ave (SE)
			& Total (\textperthousand) &  Type \uppercase\expandafter{\romannumeral1} (\%)  & Type \uppercase\expandafter{\romannumeral2} (\%) \\
			\hline
			{Lasso}&0.473(0.529); 0.720(0.510)&0.483(0.539); 0.735(0.521)& 0(0) \\
			{Dantzig} &0.501(0.533); 0.749(0.523)&0.511(0.544); 0.764(0.534)& 0(0) \\
			{MU} &0.522(0.532); 0.741(0.501)&0.532(0.543); 0.756(0.511)& 0(0)\\
			\hline
		\end{tabular}
	\end{subtable}

	\label{tab: ar1_k100}
\end{table}

\vspace{6mm}

\begin{table} [H]
	\caption{AR(4) model with $K = 100$}
	\begin{subtable}{1\textwidth}
		\sisetup{table-format=-1.2}
		\centering
		\caption{$n = 100$}
		\begin{tabular}{c|c c c }
			\hline
			Ave (SE)
			& Total (\%) &  Type \uppercase\expandafter{\romannumeral1} (\%)  & Type \uppercase\expandafter{\romannumeral2} (\%) \\
			\hline
			{Lasso}&8.18(0.45); 8.26(0.37)&2.37(0.59); 2.27(0.44)&76.0(2.13); 78.3(1.93)\\
			{Dantzig} &8.21(0.39); 8.21(0.33)&2.38(0.51); 2.21(0.43)&76.3(1.96); 78.3(2.23)
			\\
			{MU} &8.33(0.44); 8.28(0.37)&2.57(0.52); 2.33(0.44)&75.5(1.71); 77.7(1.82)\\
			\hline
		\end{tabular}
	\end{subtable}

	\bigskip

	\begin{subtable}{1\textwidth}
		\sisetup{table-format=-1.2}
		\centering
		\caption{$n = 500$}
		\begin{tabular}{c|c c c }
			\hline
			Ave (SE)
			& Total (\%) &  Type \uppercase\expandafter{\romannumeral1} (\%)  & Type \uppercase\expandafter{\romannumeral2} (\%) \\
			\hline
			{Lasso}&4.69(0.28); 4.76(0.27)&1.22(0.36); 1.30(0.37)&45.2(3.91); 45.2(3.36)\\
			{Dantzig} &4.77(0.28); 4.81(0.27)&1.21(0.45); 1.30(0.38)&46.4(3.87); 45.9(3.61)
			\\
			{MU} &4.74(0.24); 4.79(0.25)&1.18(0.37); 1.30(0.38)&46.3(3.77); 45.6(3.76)
			\\
			\hline
		\end{tabular}
	\end{subtable}

	\bigskip

	\begin{subtable}{1\textwidth}
		\sisetup{table-format=-1.2}
		\centering
		\caption{$n = 1000$}
		\begin{tabular}{c|c c c }
			\hline
			Ave (SE)
			& Total (\%) &  Type \uppercase\expandafter{\romannumeral1} (\%)  & Type \uppercase\expandafter{\romannumeral2} (\%) \\
			\hline
			{Lasso}&2.82(0.31); 2.77(0.29)&1.37(0.37); 1.34(0.35)&19.8(1.97); 19.5(2.01)\\
			{Dantzig} &2.86(0.25); 2.80(0.29)&1.35(0.30); 1.34(0.35)&20.5(1.96); 19.8(2.07)
			\\
			{MU} &2.87(0.28); 2.79(0.28)&1.37(0.35); 1.31(0.32)&20.4(2.02); 20.1(2.02)
			\\
			\hline
		\end{tabular}
	\end{subtable}

	\label{tab: ar4_k100}
\end{table}

\begin{table} [H]
	\caption{The random precision matrix model with $\alpha = 0.1$ and $K = 100$}
	\begin{subtable}{1\textwidth}
		\sisetup{table-format=-1.2}
		\centering
		\caption{$n = 100$}
		\begin{tabular}{c|c c c }
			\hline
			Ave (SE)
			& Total (\%) &  Type \uppercase\expandafter{\romannumeral1} (\%)  & Type \uppercase\expandafter{\romannumeral2} (\%) \\
			\hline
			{Lasso} &10.3(0.68); 9.84(0.65)&4.40(0.59); 3.91(0.61)& 63.9(4.19); 64.3(5.11)\\
			{Dantzig} &10.3(0.70); 9.97(0.68)&4.23(0.74); 4.08(0.56) &66.0(4.32); 64.0(4.40)\\
			{MU}&10.1(0.65); 9.96(0.69)&4.09(0.58); 4.05(0.57)&65.6(4.44); 64.0(4.24)\\
			\hline
		\end{tabular}
	\end{subtable}

	\bigskip

	\begin{subtable}{1\textwidth}
		\sisetup{table-format=-1.2}
		\centering
		\caption{$n = 500$}
		\begin{tabular}{c|c c c }
			\hline
			Ave (SE)
			& Total (\%) &  Type \uppercase\expandafter{\romannumeral1} (\%)  & Type \uppercase\expandafter{\romannumeral2} (\%) \\
			\hline
			{Lasso} &3.47(0.51); 3.65(0.50)&2.71(0.50); 2.92(0.52)&10.4(2.80); 10.2(2.90)\\
			{Dantzig} &4.02(0.51); 4.12(0.55)&3.11(0.50); 3.43(0.58)&12.2(3.06); 10.2(2.72)\\
			{MU}&4.04(0.45); 4.25(0.61)&3.12(0.50); 3.53(0.64)&12.2(3.27); 10.5(2.39)\\
			\hline
		\end{tabular}
	\end{subtable}

	\bigskip

	\begin{subtable}{1\textwidth}
		\sisetup{table-format=-1.2}
		\centering
		\caption{$n = 1000$}
		\begin{tabular}{c|c c c }
			\hline
			Ave (SE)
			& Total (\%) &  Type \uppercase\expandafter{\romannumeral1} (\%)  & Type \uppercase\expandafter{\romannumeral2} (\%) \\
			\hline
			{Lasso} &1.34(0.38); 1.47(0.32)&1.35(0.43); 1.50(0.36)&1.27(0.67); 1.32(0.71)\\
			{Dantzig} &1.77(0.38); 1.74(0.38)&1.79(0.43); 1.77(0.42)&1.74(0.87); 1.42(0.70)\\
			{MU}&1.91(0.36); 2.29(0.58)&1.92(0.40); 2.32(0.65)&1.78(0.81); 1.79(0.80)\\
			\hline
		\end{tabular}
	\end{subtable}

	\label{tab: rd0.1_k100}
\end{table}

\vspace{6mm}

\begin{table} [H]
	\caption{The random precision matrix model with $\alpha = 0.5$ and $K = 100$}
	\begin{subtable}{1\textwidth}
		\sisetup{table-format=-1.2}
		\centering
		\caption{$n = 100$}
		\begin{tabular}{c|c c c }
			\hline
			Ave (SE)
			& Total (\%) &  Type \uppercase\expandafter{\romannumeral1} (\%)  & Type \uppercase\expandafter{\romannumeral2} (\%) \\
			\hline
			{Lasso} &49.0(0.88); 49.2(0.87)&7.04(1.16); 6.08(1.04)& 90.9(1.40); 92.1(1.30)\\
			{Dantzig} &49.2(0.88); 49.2(0.87) & 6.48(0.92); 6.22(0.90) &91.9(1.10); 92.1(1.09)\\
			{MU}&49.2(0.85); 49.2(0.89) & 6.39(0.88); 6.20(0.86) & 91.9(1.11); 92.1(1.13)\\
			\hline
		\end{tabular}
	\end{subtable}

	\bigskip

	\begin{subtable}{1\textwidth}
		\sisetup{table-format=-1.2}
		\centering
		\caption{$n = 500$}
		\begin{tabular}{c|c c c }
			\hline
			Ave (SE)
			& Total (\%) &  Type \uppercase\expandafter{\romannumeral1} (\%)  & Type \uppercase\expandafter{\romannumeral2} (\%) \\
			\hline
			{Lasso} &42.5(1.30); 42.6(1.34)&12.6(1.44); 12.8(1.66)&72.2(2.64); 72.3(2.67)\\
			{Dantzig} &44.9(1.32); 44.2(1.25)&13.7(1.73); 14.4(1.72)&76.1(2.81); 73.8(2.51)\\
			{MU}&45.7(1.23); 44.7(1.29)&13.6(1.53); 14.1(1.52)&77.6(2.13); 75.2(2.21)\\
			\hline
		\end{tabular}
	\end{subtable}

	\bigskip

	\begin{subtable}{1\textwidth}
		\sisetup{table-format=-1.2}
		\centering
		\caption{$n = 1000$}
		\begin{tabular}{c|c c c }
			\hline
			Ave (SE)
			& Total (\%) &  Type \uppercase\expandafter{\romannumeral1} (\%)  & Type \uppercase\expandafter{\romannumeral2} (\%) \\
			\hline
			{Lasso} &33.0(1.43); 32.9(1.43)&13.5(1.58); 13.8(1.58)&52.6(3.17); 52.0(3.00)\\
			{Dantzig} &36.4(1.29); 35.4(1.29)&16.2(1.47); 15.8(1.82)&56.6(2.93); 55.0(3.24)\\
			{MU}&41.2(1.12); 39.7(1.26)&17.5(2.05); 17.8(1.69)&64.8(2.28); 61.7(2.25)\\
			\hline
		\end{tabular}
	\end{subtable}
	\label{tab: rd0.5_k100}
\end{table}

\section*{Acknowledgment}
This research has partially been supported by the National Science Foundation under Grant No. DMS 1713108.

\section{Appendix: proofs} \label{sec: pfs}
Recall the notation introduced in Section~\ref{subsec: uncertainty},
\begin{align*}
	&\widehat{\bm{H}} =\bm{H}+\bm{R} \\
	&\bm{\eta}^{(t)}=\bm{X}\bm{\beta}+\bm{\epsilon}^{(t)}\;\text{ for }t=1,\cdots,n \\
	&\bm{\epsilon}^{(t)}\sim N(\bm{0},\bm{\Sigma}) \;\text{ i.i.d. for }t=1,\cdots,n.
\end{align*}
In the proofs we denote by ``$c$" a positive constant that can be different in each formula.

\subsection{Proof of Lemma~\ref{lemma: net_etaDiff}}
We first consider the case $\bm{\beta}=\bm{0}$ and show the following result:
\begin{lemma} \label{lemma: net_b0}
	Let $\sigma=\max_{k\in\wb{K}}\sqrt{\sigma_{kk}}$. Under the latent block model with $\bm{\beta}=\bm{0}$, for $L>0$, if $m_{\min}\ge 16M^2\log (nK)e^{2L}$,
	%
	\begin{multline*}
		P\Big(\max_{1 \le k \le K, 1 \le t \le n} |\widehat \eta_k^{(t)} - \eta_k^{(t)}| < 8Me^L\sqrt{\log(nK)/m_k}\Big) \\
		\geq
		1 -\sqrt{\frac{2}{\pi}}\frac{nK\sigma}{\min\{L, T\}}\exp\Big(-\frac{(\min\{L, T\})^2}{2\sigma^2}\Big) - \Big(\frac{1}{nK}\Big)^{2M^2 - 1}.
	\end{multline*}
	%
\end{lemma}
\begin{proof}
	From Hoeffding's inequality we have, for any $M > 0$, that
	%
	\begin{align*}
		P\big(\sqrt{m_k} | \widetilde{p}_{k,k}^{(t)} - p_{k,k}^{(t)}| > M\,\big|\,p_{k,k}^{(t)}\big)
		& =
		P\Big(\big|\sum\nolimits_{z[i] =k, z[j] = k} (Y^{(t)}_{ij} - p_{k,k}^{(t)})/m_k\big| > {M}/{\sqrt{m_k}}\, \Big|\, p_{k,k}^{(t)}\Big)\\
		& \leq
		2\exp\big( - 2m_k \big({M}/{\sqrt{m_k}}\big)^2 \big) =
		2\exp\big(- 2M^2\big).
	\end{align*}
	%
	Thus, using the fact that, given $\bm{p} = (p_{k,k},\ldots,p_{K,K})^T$, all the $Y_{ij}$ are independent, we obtain
	%
	\begin{align*}
		P\Big( \max_{1 \leq k \leq K, 1 \leq t \leq n}\sqrt{m_k} & |\widetilde{p}_{k,k}^{(t)} - p_{k,k}^{(t)}| > M\sqrt{\log(nK)}\,\big|\,\big\{\bm{p}^{(t)}: 1\leq t\leq n\big\}\Big) \\
		& \leq
		nK\max _{1 \leq k \leq K, 1\le t \leq n}P\left( \sqrt{m_k}|\widetilde{p}_{k,k}^{(t)} - p_{k,k}^{(t)}|  > M \sqrt{\log (nK)}\, \Big|\, p_{k,k}^{(t)}\right) \\
		& \leq
		nK\exp\left(-2M^2\log(nK)\right)  =
		(nK)^{1-2M^2},
	\end{align*}
	%
	By integrating over $\{\bm{p}^{(t)}: 1\leq t\leq n\}$, we obtain
	%
	\begin{align}\label{Pineq: net_pDiff}
		P\Big( \max_{1 \leq k \leq K, 1 \leq t \leq n}\sqrt{m_k} & |\widetilde{p}_{k,k}^{(t)} - p_{k,k}^{(t)}| > M\sqrt{\log(nK)}\Big) \leq (nK)^{1-2M^2}.
	\end{align}
	%
	Note that if $|\eta_k^{(t)}| < T$, we have
	$|\widehat{\eta}_k^{(t)} - \eta_k^{(t)} |\leq |\widetilde{\eta}_k^{(t)} - \eta_k^{(t)} |$,  and we can write
	%
	\begin{align*}
		\widetilde{\eta}_k^{(t)} - \eta_k^{(t)}
		=
		\log\tfrac{\widetilde{p}^{(t)}_{k,k}}{ 1 - \widetilde{p}_{k,k}^{(t)}} - \log\tfrac{p_{k,k}^{(t)}}{1 - p_{k,k}^{(t)}}
		=
		\tfrac{1}{\xi_k^{(t)} (1-\xi_k^{(t)})}(\widetilde{p}_{k,k}^{(t)} - p_{k,k}^{(t)}),
	\end{align*}
	%
	where $\xi_{k}^{(t)}$ lies between $\widetilde{p}_{k,k}^{(t)}$ and $p_{k,k}^{(t)}$. Since $|\xi_{k}^{(t)} - p_{k,k}^{(t)}| \leq |\widetilde{p}_{k,k}^{(t)} - p_{k,k}^{(t)}|$, inequality \eqref{Pineq: net_pDiff} also applies to $|\xi_{k}^{(t)} - p_{k,k}^{(t)}|$, so that
	%
	\begin{align*}
		P\Big(\max _{1 \leq k \leq K, 1\leq t \leq n}\sqrt{m_k} |\xi_{k}^{(t)} - p_{k,k}^{(t)}| > M\sqrt{\log (nK) }\Big)
		\leq (nK)^{1-2M^2}.
	\end{align*}
	%
	It follows that
	%
	\begin{multline} \label{Pineq: net_xi}
		P\left( \xi^{(t)}_{k} <  \epsilon \text{ or }\; \xi^{(t)}_{k} > 1- \epsilon,\text{ for all }1 \leq k \leq K, 1 \leq t \leq n\right) \\
		\leq
		P\Big( p_{k,k}^{(t)} < \epsilon + \tfrac{M\sqrt{ \log (nK)}}{\sqrt{m_k}}   \text{ or }\;p_{k,k}^{(t)}  \geq 1 - \epsilon - \tfrac{M\sqrt{ \log (nK)}}{\sqrt{m_k}},\text{ for all } 1 \leq k \leq K,\;1 \leq t \leq n\Big) \\
		+ P\Big(\max _{1 \le k \le K, 1\le t \le n}\sqrt{m_k} |\xi_{k}^{(t)} - p_{k,k}^{(t)}| > M\sqrt{\log(nK)}\Big).
	\end{multline}
	%
	As for $\eta_{k}$, since $\eta_k^{(t)}-\mu_k\sim N(0,\sigma_{kk})$ and
	%
	\begin{align*}
	\max _{1 \leq k \leq K, 1\leq t \leq n}|\eta_k^{(t)}| \leq \max _{1 \leq k \leq K, 1\leq t \leq n}|\eta_k^{(t)}-\mu_k|+\max_{1\leq k\leq K}|\mu_k|,
	\end{align*}
	%
    we have, for $C > 1$, that
	%
	\begin{align*}
		P\Big(\max _{1 \leq k \leq K, 1\leq t \leq n}|\eta_k^{(t)}| > \sqrt{2\log C}+\mu_B\Big)
		& \leq
		P\Big(\max _{1 \leq k \leq K, 1\leq t \leq n}\big|(\eta_{k}^{(t)} - \mu_k)/\sqrt{\sigma_{kk}}\big| > \sqrt{2\log C}/\sigma\Big) \\
		& \leq 2nK\widetilde{\Phi}\big({\sqrt{2\log C}}/\sigma\big),
	\end{align*}
	%
	where $\widetilde{\Phi} = 1-\Phi$ ($\Phi$ is the c.d.f. of $N(0,1)$). Note that
	%
	\begin{multline*}
		P\Big(\max _{1 \leq k \leq K, 1\leq t \leq n}|\eta_k^{(t)}| > \sqrt{2\log B}\Big) \\
		=
		P\Big(\max _{1 \leq k \leq K, 1\leq t \leq n} p_{k,k}^{(t)} > \frac{1}{1 + e^{-\sqrt{2\log B}}}\Big) + P\Big(\min _{1 \leq k \leq K, 1\leq t \leq n} p_{k,k}^{(t)} < \frac{1}{1 + e^{\sqrt{2\log B}}}\Big). \quad
	\end{multline*}
	%
	Now we are using the following:
	%
	\begin{fact}
		For each $c_0>0,$ we can find $x_0=e^{1/c_0^2} > 1$ such that, for $x \ge x_0 > 1,$ we have $e^{-\sqrt{\log x}} \ge x^{-c_0}$.
	\end{fact}
	%
	\noindent
	Using this fact, it follows that, for any $L > 0,$ we have that  $\frac{1}{1 + B^{-2/L}} \geq  \frac{1}{1 + e^{-\sqrt{2\log B\,}}}$ for $B\geq e^{L^2/2}$. W.l.o.g. assume that $B > 1.$ Then,
	%
	\begin{align*}
		P\Big(\max _{1 \leq k \leq K, 1\leq t \le n} p_{k,k}^{(t)} > \frac{1}{1 + e^{-\sqrt{2\log B\,}}}\Big)
		& \geq
		P\Big(\max _{1 \leq k \leq K, 1\leq t \leq n} p_{k,k}^{(t)} > \frac{1}{1 + B^{-2/L}}\Big) \\
		& \geq
		P\Big(\max _{1 \leq k \leq K, 1\leq t \leq n} p_{k,k}^{(t)} > 1 - 2^{-1}B^{-\frac{2}{L}}\Big).
	\end{align*}
	%
	Further, using that $\frac{e^x}{1 + e^x} \ge \frac{1}{2}$ for $x > 0$, we obtain (by using the above fact again) that for any $L > 0$ and $B\geq e^{L^2/2}$,
	%
	\begin{align*}
		P\Big(\min_{1 \leq k \leq K, 1\leq t \leq n} p_{k,k}^{(t)} < \frac{1}{1 + e^{\sqrt{2\log B}}}\Big)
		& =
		P\Big(\min_{1 \leq k \leq K, 1\leq t \leq n} p_{k,k}^{(t)} < e^{-\sqrt{2\log B}} \frac{e^{\sqrt{2\log B\,}}}{1 + e^{\sqrt{2\log B\,}}}\Big)\\
		& \geq
		P\Big(\min_{1 \leq k \leq K, 1\leq t \leq n} p_{k,k}^{(t)} < \tfrac{1}{2}B^{-\frac{2}{L}}\Big).
	\end{align*}
	%
	This means that in \eqref{Pineq: net_xi} we can choose $\epsilon = \epsilon(L) = \frac{1}{4} B^{-\frac{2}{L}}$. It follows that  with this choice of $\epsilon$ (for arbitrarily large $L$) and assuming that
	%
	\begin{align}\label{ineq: net_mk}
		\max_{1 \leq k \leq K} \tfrac{M \sqrt{\log (nK)}}{\sqrt{m_k}} \leq \epsilon,
	\end{align}
	%
	we have
	%
	\begin{align*}
		P\big( \xi^{(t)}_{k} <  \epsilon \text{ or } \xi^{(t)}_{k} > 1- \epsilon,\;1 \leq k \leq K, 1 \leq t \leq n\big)
		\leq
		2nK\widetilde{\Phi}\big({\sqrt{2\log B}}/{\sigma}\big)+ (nK)^{1 - 2M^2}.
	\end{align*}
	%
	Finally, this leads to: let $\epsilon = \epsilon(L) = \frac{1}{4} B^{-\frac{2}{L}}$. If \eqref{ineq: net_mk} holds, then for $B\geq e^{L^2/2}$
	%
	\begin{align*}
		P\Big(\tfrac{1}{\xi^{(t)}_{k} (1-\xi^{(t)}_{k})} \geq {2}/{\epsilon},\;1 \leq t \leq n, 1 \leq k \leq K  \Big)
		\leq
		2nK\widetilde{\Phi}\big({\sqrt{2\log B}}/{\sigma}\big) + (nK)^{1 - 2M^2}.
	\end{align*}
	%
	Then
	%
	\begin{multline*}
		\quad P\Big(\max_{1 \leq k \leq K, 1 \leq t \leq n}\sqrt{m_k} |\widehat \eta_k^{(t)} - \eta_k^{(t)}| < 2M\epsilon^{-1}\sqrt{\log(nK)}\Big) \\
		\geq
		1 - 2nK\widetilde{\Phi}\big(\sigma^{-1}\min\big\{\sqrt{2\log B}, T\big\}\big) - ({nK})^{1 - 2M^2}, \quad
	\end{multline*}
	%
	i.e. for $L>0$, if $m_{\min}\ge 16M^2\log( nK)B^{4/L}$ and $B\geq e^{L^2/2}$,
	%
	\begin{multline*}
		\quad P\Big(\max_{1 \leq k \leq K, 1 \leq t \leq n}\sqrt{m_k} |\widehat \eta_k^{(t)} - \eta_k^{(t)}| <8MB^{\frac{2}{L}}\sqrt{\log(nK)}\Big) \\
		\geq
		1 - 2nK\widetilde{\Phi}\big(\sigma^{-1}\min\big\{\sqrt{2\log B}, T\big\}\big) - ({nK})^{1 - 2M^2}. \quad
	\end{multline*}
	%
	Choosing $B=e^{L^2/2}$, then for $L>0$,  if $m_{\min}\geq 16M^2\log (nK)e^{2L}$,
	%
	\begin{align*}
		& P\Big(\max_{1 \leq k \leq K, 1 \leq t \leq n} |\widehat \eta_k^{(t)} - \eta_k^{(t)}| <8Me^L\sqrt{{\log(nK)}/{m_{\min}}}\,\Big) \\
		& \hspace{2.2cm} \geq
		1 -\sqrt{\frac{2}{\pi}}\frac{nK\sigma}{\min\{L, T\}}\exp\Big(-\frac{(\min\{L, T\})^2}{2\sigma^2}\Big) - \Big(\frac{1}{nK}\Big)^{2M^2 - 1}.
	\end{align*}
\end{proof}
\noindent
The proof of Lemma~\ref{lemma: net_etaDiff} with $\bm{\beta} \neq \bm{0}$ is similar, but it uses
\begin{align*}
	P\big(\max _{1 \leq k \leq K, 1\leq t \leq n}|\eta_k^{(t)}| > \sqrt{2\log C}+\mu_B\big)
	\leq
	nKC^{-1/\sigma^2}/{\sqrt{\pi\log C^{1/\sigma^2}}},
\end{align*}
which comes from $\eta_k^{(t)}-\mu_k\sim N(0,\sigma_{kk})$ and
\begin{align*}
	\max _{1 \leq k \leq K, 1\leq t \leq n}|\eta_k^{(t)}| \leq \max _{1 \leq k \leq K, 1\leq t \leq n}|\eta_k^{(t)}-\mu_k|+\max_{1\leq k\leq K}|\mu_k|.
\end{align*}

\subsection{Proof of Theorem~\ref{corollary: net_mb0_asy}}
We first introduce assumption
\begin{itemize}
	\item[C0] $|\bm{R}|_\infty < \delta$ for fixed $\delta > 0$.
\end{itemize}
This assumption makes the proofs more transparent. Intermediate results are using this assumptions. When applying these intermediate results to prove the main results (that are not using assumption C0 explicitly), we will show that C0 holds with sufficiently large probability. The following fact immediately follows from the definition of $\bm{R}$ (see Section~\ref{subsec: uncertainty}):
\begin{fact} \label{claim: mb1}
	Under assumption C0, $\|\widehat{\bm{\eta}}_k-\bm{\eta}_k\|_2\leq \sqrt{n}\delta$.
\end{fact}
\noindent
We prove a series of results: Theorem~\ref{thm: lasso_type1} and \ref{thm: lasso_type2}, Corollary~\ref{corollary: mb0_asy}, which will then imply Theorem~\ref{corollary: net_mb0_asy}.  The assertion of Theorem~\ref{corollary: net_mb0_asy} then follows from Corollary~\ref{corollary: mb0_asy} together with Lemma~\ref{lemma: net_etaDiff}. The proof is an adaptation of \cite{meinshausen2006high} (proof of Theorem 1), to our more complex situation. Both of the proofs are mainly established with the property of chi-square distribution and the Lasso.  For any $\mathcal{A} \subset \wb{K}$, let the Lasso estimate $\widehat{\bm{\theta}}^{a,\mathcal{A},\lambda, lasso}$ of $\bm{\theta}^{a,\mathcal{A}}$ be defined as
\begin{align} \label{def: thetaALasso}
	\widehat{\bm{\theta}}^{a,\mathcal{A},\lambda,lasso}=\mathop{\arg\min}_{\bm{\theta}: \theta_k=0,\forall k\notin\mathcal{A}}
	\Big(n^{-1}\|(\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n) - (\widehat{\bm{H}} - \bm{1}_n\bm{\mu})\bm{\theta}\|_2^2+\lambda\|\bm{\theta}\|_1\Big).
\end{align}
\begin{claim} \label{claim: mb2}
	For problem \eqref{def: thetaALasso}, under assumption C0, for any $q>1$,
	%
	\begin{align*}
		P\Big(\|\widehat{\bm{\theta}}^{a,\mathcal{A},\lambda,lasso}\|_1\leq (q+\delta)^2\lambda^{-1}\Big)
		\geq
		1-\exp\Big(-\tfrac{q^2-\sqrt{2q^2-1}}{2}n\Big).
	\end{align*}
	%
\end{claim}
\begin{proof}
	The claim follows directly from the tail bounds of the $\chi^2$-distribution~(see \cite{Laurent2000adaptive}) and the inequality
	%
	\begin{align*}
		n\lambda\|\widehat{\bm{\theta}}^{a,\mathcal{A},\lambda,lasso}\|_1
		\leq
		\|\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n\|_2^2
		\leq
		\left(\|\widehat{\bm{\eta}}_a-\bm{\eta}_a\|_2+\|\bm{\eta}_a - \mu_a\bm{1}_n\|_2\right)^2.
	\end{align*}
	%
\end{proof}
\begin{lemma} \label{lemma: lasso_dual}
	Given $\bm{\theta}\in\mathbb{R}^{K}$, let $G(\bm{\theta})$ be a $K$-dimensional vector with elements
	\begin{align*}
		G_b(\bm{\theta}) = -2n^{-1}\langle(\widehat{\bm{\eta}}_a-\mu_a\bm{1}_n) - (\widehat{\bm{H}} - \bm{1}_n\bm{\mu}^T)\bm{\theta}, \widehat{\bm{\eta}}_b - \mu_b\bm{1}_n\rangle.
	\end{align*}
	A vector $\widehat{\bm{\theta}}$ with $\widehat{\theta}_k=0, \forall k\in\wb{K}\backslash\mathcal{A}$ is a solution to \eqref{def: thetaALasso} iff for all $b\in\mathcal{A}$, $G_b(\widehat{\bm{\theta}})=-{\rm sign}(\widehat{\theta}_b)\lambda$ in case $\widehat{\theta}_b\neq0$, and $|G_b(\widehat{\bm{\theta}})|\leq\lambda$ in case $\widehat{\theta}_b=0$. Moreover, if the solution is not unique, and $|G_b(\widehat{\bm{\theta}})|<\lambda$ for some solution $\widehat{\bm{\theta}}$, then $\widehat{\theta}_b=0$ for all solutions of \eqref{def: thetaALasso}.
\end{lemma}
\noindent
This Lemma is almost the same as Lemma A.1 in \cite{meinshausen2006high} but without normality assumption of $\widehat{\bm{H}}$. Since the Gaussian assumption is not needed, the proof is a straightforward adaptation of the proof of Lemma A.1 in \cite{meinshausen2006high}.
\begin{lemma} \label{lemma: lasso_sign}
	For every $a\in\wb{K}$, let $\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso}$ be defined as in \eqref{def: thetaALasso}. Let the penalty parameter satisfy $\lambda_n\sim dn^{-(1-\epsilon)/2}$ with some $d>0$ and $\kappa<\epsilon<\xi$. Suppose that assumptions A1 and C0 hold with  $\delta = o(n^{-(4-\xi-3\epsilon)/2})$.  Then there exists $c>0$ so that, for all $a\in\wb{K}$,
	%
	\begin{align*}
		P\big({\rm sign}(\widehat{\theta}_b^{a,{\rm ne}_a,\lambda,lasso})={\rm sign}(\theta_b^a),\; \forall b\in{\rm ne}_a\big)
		=
		1 - O\big(\exp(-cn^\epsilon)\big)
		\text{ as }n\rightarrow\infty.
	\end{align*}
	%
\end{lemma}
\begin{proof}
	Using similar notation as in \cite{meinshausen2006high}, we set
	%
	\begin{align} \label{def: thetaNeLasso}
		\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso}=\mathop{\arg\min}\limits_{\bm{\theta}: \theta_k=0,\forall k\notin{\rm ne}_a}\Big(n^{-1}\|(\widehat{\bm{\eta}}_a-\mu_a\bm{1}_n) - (\widehat{\bm{H}} - \bm{1}_n{\bm{\mu}}^T)\bm{\theta}\|_2^2+\lambda\|\bm{\theta}\|_1\Big),
	\end{align}
	%
	and for all $a,b\in\wb{K}$ with $b\in{\rm ne}_a$, we let
	%
	\begin{align} \label{def: thetaOmegaLasso}
		\widetilde{\bm{\theta}}^{a,b,\lambda}(\omega)=\mathop{\arg\min}\limits_{\bm{\theta}\in\Theta_{a,b}(\omega)}\Big(n^{-1}\|(\widehat{\bm{\eta}}_a-\mu_a\bm{1}_n) - (\widehat{\bm{H}} - \bm{1}_n{\bm{\mu}}^T)\bm{\theta}\|_2^2+\lambda\|\bm{\theta}\|_1\Big),
	\end{align}
	%
	where
	%
	\begin{align*}
		\Theta_{a,b}(\omega)=\big\{\bm{\theta}\in\mathbb{R}^{K(n)}:\theta_b=\omega;\theta_k=0,\forall k\notin{\rm ne}_a\big\}.
	\end{align*}
	%
	Setting $\omega=\widehat{\theta}_b^{a,{\rm ne}_a,\lambda,lasso}$, then  $\widetilde{\bm{\theta}}^{a,b,\lambda}(\omega)= \widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso}$,  and by Claim~\ref{claim: mb2} with $q=2$,
	%
	\begin{align*}
		P\Big(|\widehat{\theta}_b^{a,{\rm ne}_a,\lambda,lasso}|\leq(2+\delta)^2\lambda^{-1}\Big)
		\geq
		1-\exp\left(-2^{-1}n\right).
	\end{align*}
	%
	Thus, if ${\rm sign}(\widehat{\theta}_b^{a,{\rm ne}_a,\lambda,lasso})\neq{\rm sign}(\theta_b^a)$,  with probability at least $1-\exp(2^{-1}n)$, there would exist some $\omega$ with $|\omega| \leq (2+\delta)^2\lambda^{-1}$ so that  $\widetilde{\bm{\theta}}^{a,b,\lambda}(\omega)$ is a solution to \eqref{def: thetaNeLasso} but ${\rm sign}(\omega){\rm sign}(\theta_b^a)\leq0$.
	Without loss of generality, we assume $\theta_b^a>0$ since ${\rm sign}(\theta_b^a)\neq0$ for all $b\in{\rm ne}_a$. Note that by Lemma~\ref{lemma: lasso_dual}, $\widetilde{\bm{\theta}}^{a,b,\lambda}(\omega)$ can be a solution to \eqref{def: thetaNeLasso} only if $G_b(\widetilde{\bm{\theta}}^{a,b,\lambda}(\omega))\geq-\lambda$ when $\omega\leq0$. This means
	%
	\begin{multline} \label{Pineq: sign}
		\quad P\Big({\rm sign}(\widehat{\theta}_b^{a,{\rm ne}_a,\lambda,lasso})\neq{\rm sign}(\theta_b^a)\Big) \\
		\leq P\Big(\sup_{-\lambda^{-1}(2+\delta)^2\leq\omega\leq0}G_b(\widetilde{\bm{\theta}}^{a,b,\lambda}(\omega))\geq-\lambda\Big)+\exp\left(-2^{-1}n\right).
	\end{multline}
	%
	Let $\bm{r}_a^{\lambda}(\omega)=(\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n) - (\widehat{\bm{H}} - \bm{1}_n{\bm{\mu}}^T)\widetilde{\bm{\theta}}^{a,b,\lambda}(\omega)$
	and write $\widehat{\eta}_b$, $\eta_b$ as
	%
	\begin{align*}
		\widehat{\eta}_b - \mu_b = \sum_{k\in{\rm ne}_a\backslash\{b\}}\theta_k^{b,{\rm ne}_a\backslash\{b\}}(\widehat{\eta}_k - \mu_k) + \widehat{w}_b
	\end{align*}
	%
	and
	%
	\begin{align*}
		\eta_b - \mu_b = \sum_{k\in{\rm ne}_a\backslash\{b\}}\theta_k^{b,{\rm ne}_a\backslash\{b\}}(\eta_k - \mu_k) + w_b;
	\end{align*}
	%
	and $\widehat{\eta}_a$, $\eta_a$ as

	\begin{align*}
		\widehat{\eta}_a - \mu_a = \sum_{k\in{\rm ne}_a}\theta_k^a(\widehat{\eta}_k - \mu_k) + \widehat{v}_a
		\quad\text{and}\quad
		\eta_a - \mu_a = \sum_{k\in{\rm ne}_a}\theta_k^a(\eta_k - \mu_k) + v_a,
	\end{align*}
	%
	where $v_a$ and $w_b$ are independent normally distributed random variables with variances $\sigma_{-a}^2$ and $\sigma_{w,b}^2$, respectivley, and $0<\upsilon^2\leq\sigma_{-a}^2, \sigma_{w,b}^2\leq1$ by A1.2. Now we can write
	%
	\begin{align} \label{eq: etaVW}
		\eta_a - \mu_a = \sum_{k\in{\rm ne}_a\backslash\{b\}}(\theta_k^a+\theta_b^a\theta_k^{b,{\rm ne}_a\backslash\{b\}})(\eta_k - \mu_k) + \theta_b^aw_b+v_a.
	\end{align}
	%
	As in \cite{meinshausen2006high}, split the $n$-dimensional vector $\widehat{\bm{w}}_b$ of observations of $\widehat{w}_b$,  and also the vector $\bm{w}_b$ of observations of $w_b$ into the sum of two vectors, respectively,
	%
	\begin{align*}
		\widehat{\bm{w}}_b=\widehat{\bm{w}}_b^{\bot}+\widehat{\bm{w}}_b^{||}
		\quad\text{and}\quad
		\bm{w}_b=\bm{w}_b^{\bot}+\bm{w}_b^{||},
	\end{align*}
	%
	where $\bm{w}_b^{||}$ and $\widehat{\bm{w}}_b^{||}$ are contained in the at most $(|{\rm ne}_a|-1)$-dimensional space $\mathbb{W}^{||}$ spanned by the vectors $\{\bm{\eta}_k: k\in{\rm ne}_a\backslash\{b\}\}$, while  $\bm{w}_b^{\bot}$ and $\widehat{\bm{w}}_b^{\bot}$ are contained in the orthogonal complement $\mathbb{W}^{\bot}$ of $\mathbb{W}^{||}$ in $\mathbb{R}^n$. Following the proof of \cite{meinshausen2006high} (Appendix, Lemma A.2), we have
	%
	\begin{align} \label{ineq: G}
		G_b(\widetilde{\bm{\theta}}^{a,b,\lambda}(\omega))
		\leq
		-2n^{-1}\langle\bm{r}_a^{\lambda}(\omega),\widehat{\bm{w}}_b\rangle+\lambda\vartheta,
	\end{align}
	%
	where $2n^{-1}\langle\bm{r}_a^{\lambda}(\omega),\widehat{\bm{w}}_b\rangle$ can be written as $2n^{-1}\langle\bm{r}_a^{\lambda}(\omega),\widehat{\bm{w}}_b^{\bot}\rangle+2n^{-1}\langle\bm{r}_a^{\lambda}(\omega),\widehat{\bm{w}}_b^{||}\rangle$.
	By definition of $\bm{r}_a^{\lambda}(\omega)$, the orthogonality property of $\widehat{\bm{w}}_b^{\bot}$, and \eqref{eq: etaVW},
	%
	\begin{align} \label{ineq: main_Psign}
		2n^{-1}\langle\bm{r}_a^{\lambda}(\omega),\widehat{\bm{w}}_b^{\bot}\rangle
		&=
		2n^{-1}\langle(\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n) - (\widehat{\bm{H}} - \bm{1}_n{\bm{\mu}}^T)\widetilde{\bm{\theta}}^{a,b,\lambda}(\omega),\widehat{\bm{w}}_b^{\bot}\rangle\nonumber\\
		&=
		2n^{-1}\langle(\bm{\eta}_a - \mu_a\bm{1}_n)-(\bm{H} - \bm{1}_n{\bm{\mu}}^T)\widetilde{\bm{\theta}}^{a,b,\lambda}(\omega),\widehat{\bm{w}}_b^{\bot}\rangle\nonumber\\
		&\qquad\quad + 2n^{-1}\langle(\widehat{\bm{\eta}}_a-\bm{\eta}_a)-(\widehat{\bm{H}}-\bm{H})\widetilde{\bm{\theta}}^{a,b,\lambda}(\omega),\widehat{\bm{w}}_b^{\bot}\rangle\nonumber\\
		&=
		2n^{-1}(\theta_b^a-\omega)\langle\bm{w}_b^{\bot},\widehat{\bm{w}}_b^{\bot}\rangle+2n^{-1}\langle\bm{v}_a,\widehat{\bm{w}}_b^{\bot}\rangle\nonumber\\
		&\qquad\quad + 2n^{-1}\langle(\widehat{\bm{\eta}}_a-\bm{\eta}_a)-(\widehat{\bm{H}}-\bm{H})\widetilde{\bm{\theta}}^{a,b,\lambda}(\omega),\widehat{\bm{w}}_b^{\bot}\rangle\nonumber\\
		&\geq
		2n^{-1}(\theta_b^a-\omega)\langle\bm{w}_b^{\bot},\widehat{\bm{w}}_b^{\bot}\rangle-|2n^{-1}\langle\bm{v}_a,\widehat{\bm{w}}_b^{\bot}\rangle|\nonumber\\
		&\qquad\quad
		- 2n^{-1}\|(\widehat{\bm{\eta}}_a-\bm{\eta}_a)-(\widehat{\bm{H}}-\bm{H})\widetilde{\bm{\theta}}^{a,b,\lambda}(\omega)\|_2\|\widehat{\bm{w}}_b^{\bot}\|_2.
	\end{align}
	%
	\begin{claim} \label{claim: mb3}
		$\|\widehat{\bm{w}}_b-\bm{w}_b\|_2\leq (\vartheta+1)\sqrt{n}\delta$ under assumption C0.
	\end{claim}
	%
	\begin{proof}
		By assumption $\|\bm{\theta}^{b,{\rm ne}_a\backslash\{b\}}\|_1\leq\vartheta$, an application of the triangle inequality and Claim~\ref{claim: mb1} gives the assertion.
	\end{proof}

	In order to estimate the second term $|2n^{-1}\langle\bm{v}_a,\widehat{\bm{w}}_b^{\bot}\rangle|$, we first consider $|2n^{-1}\langle\bm{v}_a,\bm{w}_b^{\bot}\rangle|,$ which has already been estimated in \cite{meinshausen2006high}: for every $g>0$, there exists some $c=c(g,d)>0$ so that,
	%
	\begin{align} \label{Pineq: VW1}
		P\left(|2n^{-1}\langle\bm{v}_a,\bm{w}_b^{\bot}\rangle|\geq g\lambda\right)
		\leq
		P\left(|2n^{-1}\langle\bm{v}_a,\bm{w}_b\rangle|\geq g\lambda\right)
		= O(\exp(-cn^{\epsilon})) \quad \text{as }n\rightarrow \infty.
	\end{align}
	%
	Then for the difference $||\langle\bm{v}_a,\widehat{\bm{w}}_b^{\bot}\rangle|-|\langle\bm{v}_a,\bm{w}_b^{\bot}\rangle||$, we have
	%
	\begin{align} \label{Pineq: VW2}
		P\left(2n^{-1}| |\langle\bm{v}_a,\widehat{\bm{w}}_b^{\bot}\rangle|-|\langle\bm{v}_a,\bm{w}_b^{\bot}\rangle||\leq 4(\vartheta+1)\delta\right)
		\geq
		1 - \exp(-2^{-1}n),
	\end{align}
	%
	which follows from the inequality
	%
	\begin{align*}
		\big||\langle\bm{v}_a,\widehat{\bm{w}}_b^{\bot}\rangle|-|\langle\bm{v}_a,\bm{w}_b^{\bot}\rangle|\big|
		\leq
		\|\bm{v}_a\|_2\|\widehat{\bm{w}}_b^{\bot}-\bm{w}_b^{\bot}\|_2
		\leq
		(\vartheta+1)\sqrt{n}\delta\|\bm{v}_a\|_2
	\end{align*}
	%
	together with $\|\bm{v}_a\|_2\sim\sigma_{-a}\sqrt{\chi_n^2}$.
	Thus, by \eqref{Pineq: VW1} and \eqref{Pineq: VW2},
	%
	\begin{align} \label{Pineq: VW}
		P\left(|2n^{-1}\langle\bm{v}_a,\widehat{\bm{w}}_b^{\bot}\rangle|\geq g\lambda+4(\vartheta+1)\delta\right)
		= O(\exp(-cn^{\epsilon}))
		\quad\text{as }n\rightarrow\infty.
	\end{align}
	%
	Similarly, we have
	%
	\begin{align*}
		P\left(2n^{-1}|\langle\bm{w}_b^{\bot},\widehat{\bm{w}}_b^{\bot}\rangle-\langle\bm{w}_b^{\bot},\bm{w}_b^{\bot}\rangle|\leq 4(\vartheta+1)\delta\right)
		\geq
		1-\exp\left(-2^{-1}n\right).
	\end{align*}
	%
	Note that $\sigma_{w,b}^{-2}\langle\bm{w}_b^{\bot},\bm{w}_b^{\bot}\rangle$ follows a $\chi^2_{n-|{\rm ne}_a|+1}$ distribution for $n \ge |{\rm ne}_a|$. Using again the tail bound of the $\chi^2$-distribution from \cite{Laurent2000adaptive}, we obtain with assumption A.1.3.(a) and $\sigma_{w,b}^2\geq \upsilon^2$, that there exists $n_0$ so that for $n>n_0$,
	%
	\begin{align*}
		P(2n^{-1}\langle\bm{w}_b^{\bot},\bm{w}_b^{\bot}\rangle>\upsilon^2)
		\geq
		1-\exp\left(-32^{-1}n\right).
	\end{align*}
	%
	It follows that,
	%
	\begin{align} \label{Pineq: W}
		P\left(2n^{-1}\langle\bm{w}_b^{\bot},\widehat{\bm{w}}_b^{\bot}\rangle>\upsilon^2-4(\vartheta+1)\delta\right) = O(\exp(-cn^{\epsilon}))
		\quad\text{as }n\rightarrow\infty.
	\end{align}
	%
	For the third term of \eqref{ineq: main_Psign}, note that by definition of $\widetilde{\bm{\theta}}^{a,b,\lambda}(\omega)$,
	%
	\begin{align*}
		\|(\widehat{\bm{\eta}}_a-\bm{\eta}_a)-(\widehat{\bm{H}}-\bm{H})\widetilde{\bm{\theta}}^{a,b,\lambda}(\omega)\|_2
		&\leq
		\|\widehat{\bm{\eta}}_a-\bm{\eta}_a\|_2+\sum_{k\in{\rm ne}_a}|\widetilde{\theta}^{a,b,\lambda}_k(\omega)|\|\widehat{\bm{\eta}}_k-\bm{\eta}_k\|_2 \nonumber \\
		&\leq
		\sqrt{n}\delta\big(1+\|\widetilde{\theta}^{a,b,\lambda}(\omega)\|_1\big)
	\end{align*}
	%
	and we also have
	%
	\begin{align*}
		\|\widetilde{\bm{\theta}}^{a,b,\lambda}(\omega)\|_1-|\omega|
		& \leq
		n^{-1}\|(\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n)-\omega(\widehat{\bm{\eta}}_b - \mu_b\bm{1}_n)\|_2^2\lambda^{-1} \\
		& \leq
		\left(n^{-1/2}\|\bm{\eta}_a - \mu_a\bm{1}_n\|_2+|\omega|n^{-1/2}\|\bm{\eta}_b - \mu_b\bm{1}_n\|_2+\delta(1+|\omega|)\right)^2\lambda^{-1}.
	\end{align*}
	%
	Together with $\|\widehat{\bm{w}}_b\|_2\leq\|\bm{w}_b\|_2+\|\widehat{\bm{w}}_b-\bm{w}_b\|
	\leq
	\|\bm{w}_b\|_2+(\vartheta+1)\sqrt{n}\delta$, and the property of the $\chi^2$-distribution, we have with probability at least
	$1-3\exp\left(-2^{-1}n\right)$,
	%
	\begin{multline} \label{ineq: Weta}
		2n^{-1}\|(\widehat{\bm{\eta}}_a-\bm{\eta}_a)-(\widehat{\bm{H}}-\bm{H})\widetilde{\bm{\theta}}^{a,b,\lambda}(\omega)\|_2\|\widehat{\bm{w}}_b\|_2\\
		\leq 2\left(1+|\omega|+(1+|\omega|)^2(2+\delta)^2\lambda^{-1}\right)\left(2+(\vartheta+1)\delta\right)\delta.
	\end{multline}
	%
	Using \eqref{ineq: main_Psign}, \eqref{Pineq: VW}, \eqref{Pineq: W} and \eqref{ineq: Weta}, we obatin that with probability $1-O(\exp(-cn^\epsilon))$, as $n\rightarrow\infty$,
	%
	\begin{multline*}
		2n^{-1}\langle\bm{r}_a^{\lambda}(\omega),\widehat{\bm{w}}_b^{\bot}\rangle
		\geq
		(\theta_b^a-\omega)\left(\upsilon^2-4(\vartheta+1)\delta\right)-g\lambda-4(\vartheta+1)\delta \\
		-2\left(1+|\omega|+(1+|\omega|)^2(2+\delta)^2\lambda^{-1}\right)\left(2+(\vartheta+1)\delta\right)\delta.
	\end{multline*}
	%
	Moreover, as will be shown in Lemma~\ref{lemma: lasso_ineq}, there exists $n_g=n(g)$ so that, for all $n\geq n_g$,
	%
	\begin{multline*}
		P\Big(\inf_{\omega\leq0}\{2n^{-1}\langle\bm{r}_a^{\lambda}(\omega),\widehat{\bm{w}}_b^{||}\rangle/(1+|\omega|)\}
		\geq
		-2\left(g\lambda+(\vartheta+1)\delta\right)(2+\delta)\Big) \\
		\geq
		1-2\exp\left(-2^{-1}n\right) - \exp\left(-4^{-1}g^2n\lambda^2\right).
	\end{multline*}
	%
	Thus, with probability $1-O(\exp(-cn^\epsilon))$, as $n\rightarrow\infty$,
	%
	\begin{align} \label{ineq: RW}
		2n^{-1}\langle\bm{r}_a^{\lambda}(\omega),\widehat{\bm{w}}_b\rangle
		\geq
		& (\theta_b^a-\omega)\left(\sigma_b^2-4(\vartheta+1)\delta\right)-g\lambda-4(\vartheta+1)\delta \nonumber\\
		-&
		2\left(1+|\omega|+(1+|\omega|)^2(2+\delta)^2\lambda^{-1}\right)\left(2+(\vartheta+1)\delta\right)\delta \nonumber\\
		- & 2\left(1+|\omega|\right)\left(g\lambda+(\vartheta+1)\delta\right)(2+\delta).
	\end{align}
	%
	Note that $\lambda\sim dn^{-(1-\epsilon)/2}$ with $\epsilon < \xi$, and by A1.2 and A1.4, we have
	%
	\begin{align*}
		|\theta_b^a|
		=
		|\pi_{ab}|\sqrt{\frac{{\rm Var}(\eta_b|\bm{\eta}_{\wb{K}\backslash\{b\}})}{{\rm Var}(\eta_a|\bm{\eta}_{\wb{K}\backslash\{a\}})}} \geq v\pi_{ab} \geq cvn^{-(1-\xi)/2}.
	\end{align*}
	%
	Together with \eqref{ineq: RW}, for $\delta=o(n^{-(4-\xi-3\epsilon)/2})$, we have for any $l>0$ that
	%
	\begin{align*}
		P\Big(\inf_{-(2+\delta)^2\lambda^{-1}\leq\omega\leq0}\{2n^{-1}\langle\bm{r}_a^{\lambda}(\omega),\widehat{\bm{w}}_b\rangle\}>l\lambda\Big)
		= 1 - O(\exp(-cn^\epsilon))
		\quad\text{ as }n\rightarrow\infty.
	\end{align*}
	%
	Choosing $l=\vartheta+1$ and using \eqref{ineq: G}, we have
	%
	\begin{align*}
		P\Big(\sup_{-(2+\delta)^2\lambda^{-1}\leq\omega\leq0}G_b(\widetilde{\bm{\theta}}^{a,b,\lambda}(\omega))<-\lambda\Big)
		= 1-O(\exp(-cn^\epsilon))
		\quad \text{as }n\rightarrow\infty.
	\end{align*}
	%
	Then, by Bonferroni's inequality, assumption A1.3.(a) and \eqref{Pineq: sign},
	%
	\begin{align*}
		P\big({\rm sign}(\widehat{\theta}_b^{a,{\rm ne}_a,\lambda,lasso})={\rm sign}(\theta_b^a),\,\forall b\in{\rm ne}_a\big)
		=
		1 - O\big(\exp(-cn^\epsilon)
		\quad \text{as }n\rightarrow\infty.
	\end{align*}
	%
\end{proof}

\begin{lemma} \label{lemma: lasso_ineq}
	Under assumption C0, for any $g>0$, there exists $n_g=n(g)$ so that, for all $n\geq n_g$,
	%
	\begin{multline*}
		P\Big(\inf_{\omega\leq0}\{2n^{-1}\langle\bm{r}_a^{\lambda}(\omega),\widehat{\bm{w}}_b^{||}\rangle/(1+|\omega|)\}
		\geq
		-2\left(g\lambda+(\vartheta+1)\delta\right)(2+\delta)\Big) \\
		\geq
		1-2\exp\left(-2^{-1}n\right) - \exp\left(4^{-1}g^2n\lambda^2\right).
	\end{multline*}
	%
\end{lemma}
\begin{proof}
	Again following similar arguments as in \cite{meinshausen2006high} (Appendix, Lemma A.3), we have
	%
	\begin{align*}
		|2n^{-1}\langle\bm{r}_a^{\lambda}(\omega),\widehat{\bm{w}}_b^{||}\rangle|/(1+|\omega|)
		\leq
		2n^{-1/2}\|\widehat{\bm{w}}_b^{||}\|_2\frac{n^{-1/2}\|\bm{r}_a^{\lambda}(\omega)\|_2}{1+|\omega|}
	\end{align*}
	%
	and
	%
	\begin{align*}
		P\Big(\sup_{\omega\in\mathbb{R}}\frac{n^{-1/2}\|\bm{r}_a^{\lambda}(\omega)\|_2}{1+|\omega|} > 2+\delta\Big)
		&\leq
		P\left(n^{-1/2}\max\{\|\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n\|_2,\|\widehat{\bm{\eta}}_b - \mu_b\bm{1}_n\|_2\}>2+\delta\right)\nonumber\\
		&\leq
		2\exp\left(-2^{-1}n\right).
	\end{align*}
	%
	The last inequality uses that $\|\bm{\eta}_k - \mu_k\bm{1}_n\|_2^2\sim \chi_n^2$ and $\|\widehat{\boldsymbol \eta}_k  - {\boldsymbol \eta}_k\|_2 \leq \sqrt{n}\delta$.
	Note that $\sigma_{w,b}^{-2}\langle\bm{w}_b^{||},\bm{w}_b^{||}\rangle$ follows a $\chi^2_{|{\rm ne}_a|-1}$ distribution for large $n$ and $|{\rm ne}_a|=o(n\lambda^2)$, and thus for any $g>0$, there exists $n_g=n(g)$ so that for all $n\geq n_g$,
	%
	\begin{align*}
		P\big(n^{-1/2}\|\bm{w}_b^{||}\|_2>g\lambda\big)
		\leq
		\exp\left(-4^{-1}g^2n\lambda^2\right).
	\end{align*}
	%
	Together with Claim~\ref{claim: mb3}, for any $g>0$, there exists $n_g=n(g)$ so that, for all $n\geq n_g$,
	%
	\begin{align*}
		P\big(n^{-1/2}\|\widehat{\bm{w}}_b^{||}\|_2 > g\lambda+(\vartheta+1)\delta\big)
		\leq
		\exp\left(-4^{-1}g^2n\lambda^2\right),
	\end{align*}
	%
	and thus,
	%
	\begin{multline*}
		P\Big(\sup_{\omega\in\mathbb{R}}\big\{|2n^{-1}\langle\bm{r}_a^{\lambda}(\omega),\widehat{\bm{w}}_b^{||}\rangle|/\left(1+|\omega|\right)\big\}
		\leq
		2\left(g\lambda+(\vartheta+1)\delta\right)(2+\delta)\Big) \\
		\geq
		1 - \exp\left(-4^{-1}g^2n\lambda^2\right) - 2\exp\left(-2^{-1}n\right).
	\end{multline*}
	%
\end{proof}

\begin{theorem} \label{thm: lasso_type1}
	Assume that A1 holds and that $\bm{\mu}$ is known. Let the penalty parameter satisfy $\lambda_n\sim dn^{-(1-\epsilon)/2}$ with $d>0$ and $\kappa<\epsilon<\xi$. If, in addition, C0 holds with $\delta=o(n^{\min\{-(4-\xi-3\epsilon)/2, \epsilon-\kappa-1\}})$, then for all $a\in\wb{K}$,
	%
	\begin{align*}
		P\big(\widehat{{\rm ne}}_a^{\lambda,lasso}\subseteq{\rm ne}_a\big)
		= 1 - O\big(\exp(-cn^\epsilon)\big)
		\quad
		\text{as }n\rightarrow\infty.
	\end{align*}
	%
\end{theorem}
\begin{proof}
	Following the proof of Theorem 1, \cite{meinshausen2006high}, we have
	%
	\begin{align*}
		P\big(\widehat{{\rm ne}}_a^{\lambda,lasso}\subseteq{\rm ne}_a\big)
		= 1 - P\big(\exists b\in\wb{K}\backslash{\rm cl}_a: \widehat{\theta}_b^{a,\lambda,lasso} \neq 0\big),
	\end{align*}
	%
	and
	%
	\begin{align*}
		P\left(\exists b\in\wb{K}\backslash{\rm cl}_a:\widehat{\theta}_b^{a,\lambda,lasso}\neq0\right)
		\leq P\Big(\max_{b\in\wb{K}\backslash{\rm cl}_a}|G_b(\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso})| \geq \lambda\Big),
	\end{align*}
	%
	where
	%
	\begin{align*}
		G_b\big(\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso}\big)
		= -2n^{-1}\big\langle(\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n)-(\widehat{\bm{H}} - \bm{1}_n\bm{\mu}^T)\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso}, \widehat{\bm{\eta}}_b - \mu_b\bm{1}_n\big\rangle.
	\end{align*}
	%
	For any $b\in\wb{K}\backslash{\rm cl}_a$, write
	%
	\begin{align} \label{eq: etaNeV}
		\widehat{\eta}_b - \mu_b = \sum_{k\in{\rm ne}_a}\theta_k^{b,{\rm ne}_a}(\widehat{\eta}_k  - \mu_k) + \widehat{v}_b
		%
		\text{ and }
		%
		\eta_b - \mu_b = \sum_{k\in{\rm ne}_a}\theta_k^{b,{\rm ne}_a}(\eta_k - \mu_k) + \widetilde{v}_b,
	\end{align}
	%
	where $\widetilde{v}_b \sim N(0,\sigma_{v,b}^2)$ with $\upsilon^2\leq\sigma_{v,b}^2\leq1$ and is independent of $\{\eta_k: k\in{\rm cl}_a\}$.
	%
	\begin{claim} \label{claim: mb7}
		Under assumption C0, for any $q>1$, with probability at least $1  -(|{\rm ne}_a|+2)\exp\Big(-\frac{q^2-\sqrt{2q^2-1}}{2}n\Big)$,
		%
		\begin{multline*}
			\big|\langle\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n - (\widehat{\bm{H}} - \bm{1}_n\bm{\mu}^T)\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso}, \widehat{\bm{v}}_b\rangle - \langle\bm{\eta}_a - \mu_a\bm{1}_n - (\bm{H} - \bm{1}_n\bm{\mu}^T)\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso},\widetilde{\bm{v}}_b\rangle\big|\\
			\leq
			n\left(\lambda^{-1}(q+\delta)+1\right)\left(1+\upsilon^{-1}|{\rm ne}_a|\right)(2q+\delta)\delta;
		\end{multline*}
		%
		and with probability at least $1 -(|{\rm ne}_a|+1)\exp\Big(-\frac{q^2-\sqrt{2q^2-1}}{2}n\Big)$,
		%
		\begin{multline*}
			\big|\big(\|(\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n) - (\widehat{\bm{H}} - \bm{1}_n\bm{\mu}^T)\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso}\|_2^2-\|\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n\|_2^2\big) - \\ \big(\|(\bm{\eta}_a - \mu_a\bm{1}_n)-(\bm{H} - \bm{1}_n\bm{\mu}^T)\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso}\|_2^2-\|\bm{\eta}_a - \mu_a\bm{1}_n\|_2^2\big)\big|\\
			\leq
			n\big((\lambda^{-1}(q+\delta)+1)^2+1\big)(2q+\delta)\delta.
		\end{multline*}
		%
	\end{claim}
	\begin{proof}
		Using triangle inequality and Cauchy's inequality,
		%
		\begin{align} \label{ineq: etaDiffV}
		&
		\big|(\langle\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n) - (\widehat{\bm{H}} - \bm{1}_n\bm{\mu}^T)\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso}, \widehat{\bm{v}}_b\rangle \nonumber\\
		& \hspace{2in}
		- \langle(\bm{\eta}_a - \mu_a\bm{1}_n)- (\bm{H} - \bm{1}_n\bm{\mu}^T)\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso},\widetilde{\bm{v}}_b\rangle\big| \nonumber\\
		\leq &
		\Big(\sum_{k\in{\rm ne}_a}|\widehat{\theta}^{a,{\rm ne}_a,\lambda,lasso}_k|\|\widehat{\bm{\eta}}_k-\bm{\eta}_k\|_2+\|\widehat{\bm{\eta}}_a-\bm{\eta}_a\|_2\Big)\|\widehat{\bm{v}}_b\|_2  \nonumber\\
		& \hspace{2in} +\|\bm{\eta}_a - \mu_a\bm{1}_n) - (\bm{H} - \bm{1}_n\bm{\mu}^T)\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso}\|_2\|\widehat{\bm{v}}_b-\widetilde{\bm{v}}_b\|_2\nonumber\\
		\leq &
		\Big(\big(\sum_{k\in{\rm ne}_a}|\widehat{\theta}^{a,{\rm ne}_a,\lambda,lasso}_k|+1\big)\sqrt{n}\delta+\sum_{k\in{\rm ne}_a}|\widehat{\theta}^{a,{\rm ne}_a,\lambda,lasso}_k|\|\bm{\eta}_k - \mu_k\bm{1}_n\|_2 \nonumber\\
		& \hspace{1in} +\|\bm{\eta}_a - \mu_a\bm{1}_n\|_2\Big)\|\widehat{\bm{v}}_b-\widetilde{\bm{v}}_b\|_2
		+ \Big(\sum_{k\in{\rm ne}_a}|\widehat{\theta}^{a,{\rm ne}_a,\lambda,lasso}_k|+1\Big)\sqrt{n}\delta\|\widetilde{\bm{v}}_b\|_2.
		\end{align}
		%
		By definition of $\widehat{v}_b$ and $\widetilde{v}_b$ in \eqref{eq: etaNeV}, and by using the triangle inequality,
		%
		\begin{align*}
			\|\widehat{\bm{v}}_b-\widetilde{\bm{v}}_b\|_2
			\leq \Big(\sum_{k\in{\rm ne}_a}|\theta_k^{b,{\rm ne}_a}|+1\Big)\sqrt{n}\delta.
		\end{align*}
		%
		Moreover, by the definition of partial correlation and assumption A1.4,
		%
		\begin{align*}
			1 \geq |\pi_{bk}^{b,{\rm ne}_a}| = |\theta_k^{b,{\rm ne}_a}|\sqrt{\frac{{\rm Var}(\eta_b|\bm{\eta}_{{\rm ne}_a})}{{\rm Var}(\eta_k|\bm{\eta}_{\{b\}\cup{\rm ne}_a\backslash\{k\}})}} \geq v|\theta_k^{b,{\rm ne}_a}|,
		\end{align*}
		%
		and thus
		%
		\begin{align} \label{ineq: DiffV}
			\|\widehat{\bm{v}}_b-\widetilde{\bm{v}}_b\|_2
			\leq
			(\upsilon^{-1}|{\rm ne}_a|+1)\sqrt{n}\delta.
		\end{align}
		%
		Using \eqref{ineq: etaDiffV}, \eqref{ineq: DiffV}, Claim~\ref{claim: mb2} and the property of chi-square distribution, with probability $1  -(|{\rm ne}_a|+2)\exp\Big(-\frac{q^2-\sqrt{2q^2-1}}{2}n\Big)$,
		%
		\begin{multline*}
			\big|\langle(\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n)-(\widehat{\bm{H}} - \bm{1}_n\bm{\mu}^T)\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso}, \widehat{\bm{v}}_b\rangle - \langle(\bm{\eta}_a - \mu_a\bm{1}_n) - (\bm{H} - \bm{1}_n\bm{\mu}^T)\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso},\widetilde{\bm{v}}_b\rangle\big| \\
			\leq
			n\left((q+\delta)\lambda^{-1}+1\right)\left(1+\upsilon^{-1}|{\rm ne}_a|\right)(2q+\delta)\delta.
		\end{multline*}
		%
		The second part follows similarly.
	\end{proof}
	%
	\noindent
	Lemma~\ref{lemma: lasso_dual}, \ref{lemma: lasso_sign} and assumption A1.5 imply that, for any  $\delta=o(n^{-(4-\xi-3\epsilon)/2})$,
	there exists $c>0$ so that for all $a\in\wb{K}$ and $b\in \wb{K}\backslash{\rm cl}_a$,
	%
	\begin{multline} \label{eq: Peq_G}
		P\Big(|G_b(\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso})|
		\leq\varrho\lambda + |2n^{-1}\langle(\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n)-(\widehat{\bm{H}} - \bm{1}_n\bm{\mu}^T)\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso},\widehat{\bm{v}}_b\rangle|\Big) \\
		=
		1-O(\exp(-cn^\epsilon))
		\quad\text{as }n\rightarrow\infty.
	\end{multline}
	%
	Now we need to estimate $|2n^{-1}\langle(\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n)-(\widehat{\bm{H}} - \bm{1}_n\bm{\mu}^T)\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso},\widehat{\bm{v}}_b\rangle|$.
	Note that by Claim~\ref{claim: mb7}, for any $\delta=O(1)$, there exists some constant $B>0$ and $c>0$ so that for $n\rightarrow\infty$, with probability $1-O(\exp(-cn^\epsilon))$,
	%
	\begin{multline} \label{ineq: PetaVDiff}
		\big|\langle(\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n)-(\widehat{\bm{H}} - \bm{1}_n\bm{\mu}^T)\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso},\widehat{\bm{v}}_b\rangle - \langle(\bm{\eta}_a - \mu_a\bm{1}_n) - (\bm{H} - \bm{1}_n\bm{\mu}^T)\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso},\widetilde{\bm{v}}_b\rangle\big| \\
		\leq
		Bn^{3/2-\epsilon/2+\kappa}\delta.
	\end{multline}
	%
	We already know that $\widetilde{v}_b{\perp\!\!\!\perp}\bm{\eta}_{{\rm cl}_a}$\footnote{This follows from the Markov properties of the conditional independence graph and the contraction property of conditional independence.} and $\widehat{\bm{\eta}}_{{\rm cl}_a}{\perp\!\!\!\perp}\bm{\eta}_{\wb{K}\backslash{\rm cl}_a}|\bm{\eta}_{{\rm cl}_a}$, we have $\widetilde{v}_b{\perp\!\!\!\perp}\widehat{\bm{\eta}}_{{\rm cl}_a}|\bm{\eta}_{{\rm cl}_a}$ since $\widetilde{v}_b=\eta_b-\sum_{k\in{\rm ne}_a}\theta_k^{b,{\rm ne}_a}\eta_k$. Thus, $\widetilde{v}_b{\perp\!\!\!\perp}\{\bm{\eta}_{{\rm cl}_a}, \widehat{\bm{\eta}}_{{\rm cl}_a}\}$. Here ${\perp\!\!\!\perp}$ denotes independence. Conditional on $\{\bm{H}_{{\rm cl}_a}, \widehat{\bm{H}}_{{\rm cl}_a}\}=\{\bm{\eta}_k, \widehat{\bm{\eta}}_k: k\in{\rm cl}_a\}$, the random variable
	%
	\begin{align*}
		\langle(\bm{\eta}_a - \mu_a\bm{1}_n) - (\bm{H} - \bm{1}_n\bm{\mu}^T)\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso}, \widetilde{\bm{v}}_b\rangle
	\end{align*}
	%
	is normally distributed with mean zero and variance

	$$
	\sigma_{v,b}^2\|(\bm{\eta}_a - \mu_a\bm{1}_n) - (\bm{H} - \bm{1}_n\bm{\mu}^T)\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso}\|_2^2.
	$$
	By definition of $\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso}$,
	%
	\begin{align*}
		\|(\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n) - (\widehat{\bm{H}} - \bm{1}_n\bm{\mu}^T)\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso}\|_2
		\leq
		\|\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n\|_2;
	\end{align*}
	%
	and by Claim~\ref{claim: mb7}, for $\delta=O(1)$, there exists constant $c, B>0$ so that, with probability $1-O(\exp(-cn^\epsilon))$, as $n \rightarrow \infty$, as $n \rightarrow \infty$,
	%
	\begin{multline*}
		\big|\big(\|(\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n) - (\widehat{\bm{H}} - \bm{1}_n\bm{\mu}^T)\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso}\|_2 - \|\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n\|_2^2\big) \\
		-(\|(\bm{\eta}_a - \mu_a\bm{1}_n) - (\bm{H} - \bm{1}_n\bm{\mu}^T)\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso}\|_2^2-\|\bm{\eta}_a - \mu_a\bm{1}_n\|_2^2)\big|
		\leq
		Bn^{2-\epsilon}\delta.
	\end{multline*}
	%
	Futhermore, for $\delta=O(n^{\min\{0,2\tau+\epsilon-2\}})$, there exists $c_t>0$, such that for  $t_a = c_tn^{\tau}$,
	%
	\begin{multline*}
		P\left(\|(\bm{\eta}_a - \mu_a\bm{1}_n) - (\bm{H} - \bm{1}_n\bm{\mu}^T)\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso}\|_2^2 \leq \|\bm{\eta}_a - \mu_a\bm{1}_n\|_2^2+ t_a^2\right) \\
		= 1-O(\exp(-cn^\epsilon)) \quad \text{ as }n\rightarrow\infty,
	\end{multline*}
	%
	thus, $|2n^{-1}\langle(\bm{\eta}_a - \mu_a\bm{1}_n) - (\bm{H} - \bm{1}_n\bm{\mu}^T)\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso},\widetilde{\bm{v}}_b\rangle|$
	is stochastically smaller than $|2n^{-1}(\langle\bm{\eta}_a - \mu_a\bm{1}_n, \widetilde{\bm{v}}_b\rangle+t_az_b)|$ with probability $1-O(\exp(-cn^\epsilon))$, as $n \rightarrow \infty$, where $z_b\sim N(0,\sigma_{v,b}^2)$ and is independent of other random variables. Since $\widetilde{v}_b$ and $\eta_a$ are independent, ${\rm E}(\eta_a\widetilde{v}_b)=0$. Using the Gaussianity and Bernstein's inequality,
	%
	\begin{align*}
		P\left(|2n^{-1}(\langle\bm{\eta}_a - \mu_a\bm{1}_n, \widetilde{\bm{v}}_b\rangle+t_az_b)|\geq(1-\varrho)\lambda/2\right)
		= O\left(\exp\{-cn^{\min\{\epsilon,1+\epsilon-2\tau\}}\}\right)
		\text{ as } n\rightarrow \infty.
	\end{align*}
	%
	and thus for $\delta=O(n^{\min\{0,2\tau+\epsilon-2\}})$, as $n\rightarrow\infty$,
	%
	\begin{multline} \label{ineq: PetaV}
		P\big(|2n^{-1}\langle(\bm{\eta}_a - \mu_a\bm{1}_n) - (\bm{H} - \bm{1}_n\bm{\mu}^T)\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso},\widetilde{\bm{v}}_b\rangle|\geq(1-\varrho)\lambda/2\big) \\
		=
		O\left(\exp\{-cn^{\min\{\epsilon,1+\epsilon-2\tau\}}\}\right).
	\end{multline}
	%
	By \eqref{eq: Peq_G}, \eqref{ineq: PetaVDiff} and \eqref{ineq: PetaV}, for $\delta=o(n^{-(4-\xi-3\epsilon)/2})$ and $\delta=O(n^{2\tau+\epsilon-2})$, there exists $c, B>0$, with probability $1-O\left(\exp\{-cn^{\min\{\epsilon,1+\epsilon-2\tau\}}\}\right)$, as $n\rightarrow\infty$,
	%
	\begin{align*}
		|G_b(\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso})|<(1+\varrho)\lambda/2+Bn^{1/2-\epsilon/2+\kappa}\delta,
	\end{align*}
	%
	and we obtain that for $\delta=o(n^{\min\{-(4-\xi-3\epsilon)/2, \epsilon-\kappa-1\}})$,
	%
	\begin{align*}
		P\Big(\max_{b\in\wb{K}\backslash\{a\}}|G_b(\widehat{\bm{\theta}}^{a,{\rm ne}_a,\lambda,lasso})|\geq\lambda\Big) = O\big(\exp(-cn^\epsilon)\big)\quad\text{as }n\rightarrow\infty.\\[-30pt]
	\end{align*}
	%
\end{proof}

\begin{theorem} \label{thm: lasso_type2}
	Let assumption A1 hold and assume $\bm{\mu}$ to be known. Let the penalty parameter satisfy $\lambda_n\sim dn^{-(1-\epsilon)/2}$ with some $d>0$ and $\kappa<\epsilon<\xi$. If, in addition, C0 holds with $\delta = o(n^{\min\{-(4-\xi-3\epsilon)/2, \epsilon-\kappa-1\}})$, for all $a\in\wb{K}$,
	%
	\begin{align*}
		P({\rm ne}_a\subseteq\widehat{{\rm ne}}_a^{\lambda,lasso})=1-O(\exp(-cn^\epsilon))
		\quad
		\text{as }n\rightarrow\infty.
	\end{align*}
	%
\end{theorem}
\begin{proof}
	Using Theorem~\ref{thm: lasso_type1} and Lemma~\ref{lemma: lasso_sign}, the proof is similar to \cite{meinshausen2006high}, proof of Theorem 2.
\end{proof}

\begin{corollary} \label{corollary: mb0_asy}
	Let assumption A1 hold and assume $\bm{\mu}$ to be known. Let the penalty parameter satisfy $\lambda_n\sim dn^{-(1-\epsilon)/2}$ with some $d>0$ and $\kappa<\epsilon<\xi$. If, in addition, C0 holds with  $\delta = o(n^{\min\{-(4-\xi-3\epsilon)/2, \epsilon-\kappa-1\}})$, then there exists $ c>0$ so that
	\begin{align*}
		P(\widehat{E}^{\lambda,lasso}=E)=1- O(\exp(-cn^\epsilon))
		\quad
		\text{as }n\rightarrow\infty.
	\end{align*}
\end{corollary}
\begin{proof}
	Note that $\widehat{E}^{\lambda, lasso}\neq E$ if and only if there exists $a\in\wb{K}$ so that $\widehat{{\rm ne}}_a^{\lambda, lasso}\neq {\rm ne}_a$. The result now follows from Theorem~\ref{thm: lasso_type1} and Theorem~\ref{thm: lasso_type2} by using Bonferroni's inequality and assumption A1.1.
\end{proof}

\subsection{Proof of Theorem~\ref{corollary: net_ds0_asy}}
We prove a series of results, which will then imply Theorem~\ref{corollary: net_ds0_asy}. The asseration of Theorem~\ref{corollary: net_ds0_asy} follows from Corollary~\ref{thm: ds0_asy} together with Lemma~\ref{lemma: net_etaDiff}. First we introduce some notation. Let
\begin{align*}
	& s=\max_{a\in\wb{K}}|{\rm ne}_a| \\
	& \bm{\Psi}=\frac{1}{n}(\bm{H} - \bm{1}_n\bm{\mu}^T)^T(\bm{H} - \bm{1}_n\bm{\mu}^T) \\
	& \bm{\Psi}^a=\frac{1}{n}(\bm{H}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T)^T(\bm{H}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T) \quad\text{for each }a\in\wb{K}.
\end{align*}
We also define, for each $a\in\wb{K}$, $1\leq s \leq K$ and $1\leq q\leq \infty$
\begin{align} \label{def: kappa}
	\kappa^a_q(s)=\min_{J:|J|\leq s}\Big(\min_{\bm{\Delta}\in\mathbb{R}^p: \|\bm{\Delta}_{J^c}\|_1\leq\|\bm{\Delta}_{J}\|_1\atop \|\bm{\Delta}\|_q=1 }\big|\bm{\Psi}^a\bm{\Delta}\big|_\infty\Big).
\end{align}
\begin{claim} \label{claim: ds2}
	For $n$ large enough, $\lvert\bm{\Psi}-\bm{\Sigma}\rvert_\infty<\frac{1}{2\alpha s}$ with probability no less than $1-K^2\exp\left(-\frac{n}{150\alpha^2s^2\sigma^4}\right)$.
\end{claim}
\begin{proof}
	Denote $\bm{\Sigma}=(\sigma_{ab})_{K\times K}$.
	For $b\neq a$, we have $\bm{\eta}_b - \mu_b\bm{1}_n=\sigma_{ba}\sigma_{aa}^{-1}(\bm{\eta}_a - \mu_a\bm{1}_n)+\bm{v}_{ba}$, where each element of $\bm{v}_{ba}$ is a zero mean normal with variance $\sigma_{b|a}\leq\sigma_{bb}$; and
	%
	\begin{align*}
	& \; P\Big(\Big|\frac{1}{n}(\bm{\eta}_a - \mu_a\bm{1}_n)^T(\bm{\eta}_b - \mu_b\bm{1}_n)-\sigma_{ab}\Big|>\frac{1}{2\alpha s}\Big) \\
	\leq &\;
	P\left(\sigma_{ba}\Big|\frac{(\bm{\eta}_a - \mu_a\bm{1}_n)^T(\bm{\eta}_a - \mu_a\bm{1}_n)}{n\sigma_{aa}}-1\Big|>\frac{1}{3\alpha s}\right) + P\left(\Big|\frac{(\bm{\eta}_a - \mu_a\bm{1}_n)^T\bm{v}_{ba}}{n}\Big|>\frac{1}{6\alpha s}\right)
	\end{align*}
	%
	Using a tail bound for the chi-squared distribution from \cite{Johnstone2001chisquare}, Bernstein's inequality and Bonferroni's inequality, we have for $n$ large enough,
	%
	\begin{align*}
		P\Big(\lvert\bm{\Psi}-\bm{\Sigma}\rvert_\infty>\frac{1}{2\alpha s}\Big)
		\leq
		K^2\exp\left(-\frac{n}{150\alpha^2s^2\sigma^4}\right).
	\end{align*}
	%
\end{proof}

\begin{claim} \label{claim: ds3}
	Under assumption C0, for any $q >0$, we have with probability no less than $1-2K^2\exp\big(-\frac{nq^2}{4\sigma^4+2\sigma^2q}\big) - 2K\exp\left(-2^{-1}n\right)$,
	%
	\begin{align*}
		\max_{a\in\wb{K}}\big|n^{-1}(\widehat{\bm{H}}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T)^T\bm{\xi}_a\big|_\infty
		\leq
		q+4\sigma\delta+\delta^2.
	\end{align*}
	%
\end{claim}
\begin{proof}
	By straightforward calculation,
	%
	\begin{align*}
		&\max_{a\in\wb{K}}\big|n^{-1}(\widehat{\bm{H}}_{-a} -  \bm{1}_n\bm{\mu}_{-a}^T)^T\bm{\xi}_a\big|_\infty \\
		\leq
		&
		\max_{a\in\wb{K}}\big|n^{-1}(\widehat{\bm{H}}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T)^T\bm{v}_a\big|_\infty + \max_{a\in\wb{K}}\big|n^{-1}(\widehat{\bm{H}}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T)^T(\widehat{\bm{\eta}}_a-\bm{\eta}_a)\big|_\infty \\
		\leq
		&
		\max_{a\in\wb{K}}\max_{ b\in\wb{K}\backslash\{a\}}n^{-1}\left|(\bm{\eta}_{b} - \mu_b\bm{1}_n)^T\bm{v}_a\right| + n^{-1/2}\big(\max_{a\in\wb{K}}\|\bm{\eta}_a - \mu_a\bm{1}_n\|_2 + \max_{a\in\wb{K}}\|\bm{v}_a\|_2\big)\delta+\delta^2.
	\end{align*}
	%
	Note that $n^{-1}|(\bm{\eta} - \mu_a\bm{1}_n)^T_{b}\bm{v}_a|$ can be estimated by Bernstein's inequality, and both $\|\bm{\eta}_a - \mu_a\bm{1}_n\|_2^2$ and $\sigma_{aa}^{-1}\|\bm{v}_a\|_2^2$ are chi-square distributed. Thus, together with Bonferroni's inequality, gives
	%
	\begin{align*}
		P\Big(\max_{a\in\wb{K}}\max_{ b\in\wb{K}\backslash\{a\}}n^{-1}\left|(\bm{\eta}_{b} - \mu_b\bm{1}_n)^T\bm{v}_a\right|\geq q\Big)
		\leq
		2K^2\exp\left(-\frac{nq^2}{4\sigma^4+2\sigma^2 q}\right)
	\end{align*}
	%
	and
	%
	\begin{align*}
		P\Big(n^{-1/2}\big(\max_{a\in\wb{K}}\|\bm{\eta}_a-\mu_a\bm{1}_n\|_2+\max_{a\in\wb{K}}\|\bm{v}_a\|_2\big) > 4\sigma\Big)
		\leq
		2K\exp\left(-2^{-1}n\right).
	\end{align*}
	%
\end{proof}
\begin{claim} \label{claim: ds4}
	For any $\alpha > \|\bm{\Sigma}^{-1}\|_\infty$, there exists $n_\alpha=n(\alpha)$ so that, for all $n\geq n_\alpha$,
	%
	\begin{align*}
		P\Big(\max_{a\in\wb{K}}\kappa^a_\infty(s)\geq \|\bm{\Sigma}^{-1}\|_\infty^{-1}-\alpha^{-1}\Big)
		\leq
		1 - K^2\exp\left(-\frac{n}{150\alpha^2s^2\sigma^4}\right)
	\end{align*}
	%
\end{claim}
\begin{proof}
	The result follows from Claim~\ref{claim: ds2} and the inequality below: for all $a\in\wb{K}$,
	%
	\begin{align*}
		\kappa^a_\infty(s)
		& \geq
		\min_{J:|J|\leq s}\Big(\min_{\bm{\Delta}\in\mathbb{R}^p: \|\bm{\Delta}_{J^c}\|_1\leq\|\bm{\Delta}_{J}\|_1\atop |\bm{\Delta}|_\infty=1 }\left|\bm{\Sigma}\bm{\Delta}\right|_\infty\Big)
		- \max_{J: |J|\leq s}\Big(\max_{\bm{\Delta}\in\mathbb{R}^p: \|\bm{\Delta}_{J^c}\|_1 \leq \|\bm{\Delta}_{J}\|_1 \atop |\bm{\Delta}|_\infty=1 }\left|\left(\bm{\Psi}^a-\bm{\Sigma}\right)\bm{\Delta}\right|_\infty\Big) \\
		& \geq
		\|\bm{\Sigma}\|_\infty^{-1} - \max_{J: |J|\leq s}\max_{\bm{\Delta}\in\mathbb{R}^p: \|\bm{\Delta}_{J^c}\|_1 \leq \|\bm{\Delta}_{J}\|_1 \atop |\bm{\Delta}|_\infty=1 }\|\bm{\Delta}\|_1\left|\bm{\Psi}^a - \bm{\Sigma}\right|_\infty \\
		& \geq
		\|\bm{\Sigma}^{-1}\|_\infty^{-1} - 2s\left|\bm{\Psi} - \bm{\Sigma}\right|_\infty.
	\end{align*}
	%
\end{proof}

\begin{theorem} \label{thm: ds0_asy}
	Let assumption B1 hold and assume $\bm{\mu}$ to be known. Let $\lambda_n^{-1} = O(n^{\frac{1-\epsilon}{2}})$ with some $\epsilon>0$ be such that $\xi>\epsilon>2\kappa-2p+1$. If, in addition, C0 holds with  $\delta = O(n^{-p})$ for some $p>\kappa+(1-\xi)/2$,
	\begin{align*}
		P\big(\widehat{E}^{\lambda,ds} = E\big)
		= 1 - O\big(\exp(-cn^{\min\{\epsilon,1-2\kappa\}})\big)
		\quad
		\text{as }n\rightarrow\infty.
	\end{align*}
\end{theorem}
\begin{proof}
	Note that by Claim~\ref{claim: ds3}, the property of chi-square distribution,  and
	%
	\begin{align*}
		&\, \big|n^{-1}(\widehat{\bm{H}}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T)^T\big((\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n)-(\widehat{\bm{H}} - \bm{1}_n\bm{\mu}^T)\bm{\theta}^a\big)\big|_\infty \\
		= &\,
		\big|n^{-1}(\widehat{\bm{H}}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T)^T\big((\bm{H}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T)\bm{\theta}^a_{-a} + \bm{\xi}_a - (\widehat{\bm{H}}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T)\bm{\theta}^a_{-a}\big)\big|_\infty \\
		\leq &\,
		\big|n^{-1}(\widehat{\bm{H}}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T)^T\bm{R}_{-a}\bm{\theta}^a_{-a}\big|_\infty + \max_{a\in\wb{K}}\big|n^{-1}(\widehat{\bm{H}}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T)^T\bm{\xi}_a\big|_\infty\\
		\leq &\,
		n^{-1/2}\delta\big(\sqrt{n}\delta + \max_{b\in\wb{K}}\|\bm{\eta}_b - \mu_b\bm{1}_n\|_2\big)s\upsilon^{-1} + \max_{a\in\wb{K}}\big|n^{-1}(\widehat{\bm{H}}_{-a} - \bm{1}_n\bm{\mu}_{-a})^T\bm{\xi}_a\big|_\infty,
	\end{align*}
	%
	we have for each $q > 0$, that with probability no less than $1 - 2K\exp\left(-2^{-1}n\right) - 2K^2\exp\left(-\tfrac{nq^2}{4\sigma^4+2\sigma^2q}\right)$,
	%
	\begin{multline*}
		\max_{a\in\wb{K}}\big|n^{-1}(\widehat{\bm{H}}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T)^T\big((\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n) - (\widehat{\bm{H}} - \bm{1}_n\bm{\mu}^T)\bm{\theta}^a\big)\big|_\infty \\
		\leq q+(4\sigma+\delta)\delta+(2\sigma+\delta)\upsilon^{-1}s\delta
	\end{multline*}
	%
	Note that $\lambda_n^{-1} = O(n^{\frac{1-\epsilon}{2}})$, which means there exits $b > 0$ so that $\lambda_n \geq bn^{-\frac{1-\epsilon}{2}}$ for all $n$ large enough. Moreover, $s\delta = O(n^{\kappa-p}) = o(\lambda_{a,n}(\|\bm{\theta}^a\|_1))$.   Choosing $q=bn^{-\frac{1-\epsilon}{2}}/2$, we obtain $q+2\upsilon^{-1}s\delta \leq \lambda_{a,n}(\|\bm{\theta}^a\|_1)$ for all large $n$. Thus
	%
	\begin{multline*}
		P\big(\big|n^{-1}(\widehat{\bm{H}}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T)^T\big((\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n) - (\widehat{\bm{H}} - \bm{1}_n\bm{\mu}^T)\bm{\theta}^a\big)\big|_\infty
		\leq \lambda_{a,n}(\|\bm{\theta}^a\|_1), \,\forall a\in\wb{K} \big) \\
		= 1 - O(\exp(-cn^{\epsilon}))
		\quad\text{as }n\rightarrow\infty,
	\end{multline*}
	%
	which means that the true parameter $\bm{\theta}^a$ falls into the feasible set of problem (3.8) with probability $1-O(\exp(-cn^\epsilon))$ as $n\rightarrow\infty$.  With $\bm{\Delta}^a=\widetilde{\bm{\theta}}^{a,\lambda,ds}_{-a}-\bm{\theta}^a_{-a}$, we obtain
	%
	\begin{align*}
	\left|\bm{\Psi}^a\bm{\Delta}^a\right|_\infty
	\leq &  \big|n^{-1}(\widehat{\bm{H}}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T)^T(\bm{H}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T)\bm{\Delta}^a\big|_\infty + \big|n^{-1}\bm{R}_{-a}^T(\bm{H}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T)\bm{\Delta}^a\big|_\infty \\
	\leq  &
	\big|n^{-1}(\widehat{\bm{H}}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T)^T((\widehat{\bm{\eta}}_a - \mu_a\bm{1}_n) - (\widehat{\bm{H}}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T)\widetilde{\bm{\theta}}_{-a}^{a,\lambda,ds})\big|_\infty  \\
	& \hspace{2cm} + \big|n^{-1}(\widehat{\bm{H}}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T)^T\bm{\xi}_a\big|_\infty + \big|n^{-1}(\widehat{\bm{H}}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T)^T\bm{R}_{-a}\widetilde{\bm{\theta}}^{a,ds}_{-a}\big|_\infty \\
	& \hspace{6cm} + \big|n^{-1}\bm{R}_{-a}^T(\bm{H}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T)\bm{\Delta}^a\big|_\infty  \\
	\leq  &
	\lambda_{a,n}\big(\|\widetilde{\bm{\theta}}^{a,\lambda,ds}\|_1\big) + \big|n^{-1}(\widehat{\bm{H}}_{-a} - \bm{1}_n\bm{\mu}_{-a}^T)^T\bm{\xi}_a\big|_\infty +\delta^2\|\widetilde{\bm{\theta}}_{-a}^{a,\lambda,ds}\|_1  \\
	& \hspace{5cm} + n^{-1/2}\big(2\|\widetilde{\bm{\theta}}_{-a}^{a,\lambda,ds}\|_1+\|\bm{\theta}_{-a}^a\|_1\big)\delta\max_{b\in\wb{K}}\|\bm{\eta}_b\|_2.
	\end{align*}
	%
	Then, by Claim~\ref{claim: ds3}, the definition of $\widetilde{\bm{\theta}}^{a,\lambda,ds}_{-a}$ and the property of chi-square distribution, the following holds. For any constant $d>0$, there exists $c=c(d)>0$ so that as $n\rightarrow\infty$,
	%
	\begin{align*}
		P\big(\left|\bm{\Psi}^a\bm{\Delta}^a\right|_\infty \leq \lambda_{a,n}(\|\widetilde{\bm{\theta}}^{a,\lambda,ds}\|_1) + dn^{-\frac{1-\epsilon}{2}} +4\upsilon^{-1}s\delta, \, \forall a\in\wb{K}\big)
		=
		1-O\left(\exp(-cn^{\epsilon})\right).
	\end{align*}
	%
	By definition, we have  $\kappa_\infty(s)|\bm{\Delta}^a|_\infty\leq \left|\bm{\Psi}^a\bm{\Delta}^a\right|_\infty$. Using Claim~\ref{claim: ds4} with assumption B1.3.(b), there exists a constant $t > 0$ and some $c > 0$ so that
	%
	\begin{align*}
		P\big(\left|\bm{\Delta}^a\right|_\infty \leq t\lambda_{a,n}(\|\widetilde{\bm{\theta}}^{a,\lambda,ds}\|_1),\, \forall a\in\wb{K}\big)
		=
		1-O\left(\exp(-cn^{\min\{1-2\kappa,\epsilon\}})\right)
		\quad\text{for }n\rightarrow\infty.
	\end{align*}
	%
	The assertion of Theorem~\ref{thm: ds0_asy} follows from the fact that $|\theta_b^a|=\Omega(n^{-(1-\xi)/2})$ for all $b\in{\rm ne}_a$, $a\in\wb{K}$, $\max_{a\in\wb{K}}\lambda_{a,n}(\|\bm{\theta}^a\|_1) = o(n^{-\frac{1-\xi}{2}}(\log n)^{-1})$, and $\log n \rightarrow\infty$.
\end{proof}

\subsection{Proof for Section~\ref{subsec: extension}}
To verify Theorem~\ref{thm: net_beta}, first observe that for all $k\in\wb{K}$, under assumption C0,
\begin{align*}
	|\widebar{\widehat{\eta}}_k-\mu_k|
	&\leq
	|\widebar{\widehat{\eta}}_k-\widebar{\eta}_k|+|\widebar{\eta}_k-\mu_k| \\
	& =
	\frac{1}{n}\|\bm{1}_n^T(\widehat{\bm{\eta}}_k-\bm{\eta}_k)\|_2+|\widebar{\eta}_k-\mu_k| \\
	&\leq
	\frac{1}{\sqrt{n}}\|\widehat{\bm{\eta}}_k-\bm{\eta}_k\|_2+|\widebar{\eta}_k-\mu_k|\leq
	\delta + |\widebar{\eta}_k-\mu_k|.
\end{align*}
Recalling that $\sqrt{n}(\widebar{\eta}_k-\mu_k)/\sqrt{\sigma_{kk}}\sim N(0,1)$ and  $\sigma=\max_{k\in\wb{K}}\sqrt{\sigma_{kk}}$, we obtain
\begin{align} \label{Pineq: maxMu}
	P\Big(\max_{k\in\wb{K}}|\widebar{\eta}_k-\mu_k| \geq \delta_\mu\Big)
	\leq
	2\sum_{k\in\wb{K}}\widetilde{\Phi}\left(\frac{\sqrt{n}\delta_\mu}{\sqrt{\sigma_{kk}}}\right)
	\leq
	\sqrt{\frac{2}{\pi}}\frac{K\sigma}{\sqrt{n}\delta_\mu}\exp\left( -\frac{n\delta_\mu^2}{2\sigma^2}\right),
\end{align}
and thus,
\begin{align*}
	P\Big(\max_{k\in\wb{K}}|\widebar{\widehat{\eta}}_k-\mu_k|\geq \delta+\delta_\mu\Big)
	\leq
	\sqrt{\frac{2}{\pi}}\frac{K\sigma}{\sqrt{n}\delta_\mu}\exp\left( -\frac{n\delta_\mu^2}{2\sigma^2}\right).
\end{align*}
Using the trivial bound $\|\widebar{\widehat{\bm{\eta}}}-\bm{\mu}\|_2\leq \sqrt{K}\max_{k\in\wb{K}}|\widebar{\widehat{\eta}}_k-\mu_k|$, we have
\begin{align*}
	P\left(n^{\frac{b-\gamma}{2}}\|\widebar{\widehat{\bm{\eta}}}-\bm{\mu}\|_2 \geq n^{\frac{b-\gamma}{2}}K^{\frac{1}{2}}(\delta+\delta_\mu)\right)
	\leq
	\sqrt{\frac{2}{\pi}}\frac{K\sigma}{\sqrt{n}\delta_\mu}\exp\left( -\frac{n\delta_\mu^2}{2\sigma^2}\right)
	\quad
	\text{as }n\rightarrow\infty.
\end{align*}
Thus, for $\delta=o(n^{-\frac{b}{2}})$ with $b<1,$ and $K=O(n^{\gamma})$, we obtain that, for any fixed $\delta_0>0$,
\begin{align*}
	P\big(n^{\frac{b-\gamma}{2}}\|\widebar{\widehat{\bm{\eta}}}-\bm{\mu}\|\geq \delta_0\big) =
	O\left(\exp(-cn^{1-b})\right)
	\quad
	\text{as }n\rightarrow\infty,
\end{align*}
which, if combined with Lemma~\ref{lemma: net_etaDiff} (note that the above arguments assume C0) is the first statement of Theorem~\ref{thm: net_beta}. Recall that  $\bm{\beta} = \bm{X}^+\bm{\mu} $ and $ \widehat\bm{\beta} = \bm{X} ^+\overline{\widehat{\bm{\eta}}}$. Using the assumptions, ${\rm rank}(\bm{X})=L$ and $\sigma_{\max}(\bm{X})=O(1)$, we also obtain that, for fixed $\delta_0>0$,
\begin{align*}
	P\big(n^{\frac{b-\gamma}{2}}\|\widehat{\bm{\beta}}-\bm{\beta}\| \geq \delta_0\big) = O\left(\exp(-cn^{1-b})\right)
	\quad
	\text{as }n\rightarrow\infty.
\end{align*}
Together with Lemma~\ref{lemma: net_etaDiff}, we have the results in Theorem~\ref{thm: net_beta}.

For Corollary~\ref{corollary: net_mb_asy} and Corollary~\ref{corollary: net_ds_asy}, recall that, since $\bm{\mu}$ is unknown, we replace $\widehat{\bm{H}}$ by $\widehat{\bm{H}}-\bm{1}_n\widebar{\widehat{\bm{\eta}}}^{T}$ in \eqref{def: thetaLasso}. For our model (cf. \ref{model: matrix}), this means that we replace $\bm{R}$ by
\begin{align} \label{def: tildeR}
	\widetilde{\bm{R}}=(\widehat{\bm{H}} - \bm{1}_n\widebar{\widehat{\bm{\eta}}}^{T}) - (\bm{H} - \bm{1}_n{\bm{\mu}}^T) = \left((\widehat{\eta}_k^{(t)}-\widebar{\widehat{\eta}}_k)-(\eta_k^{(t)}-\mu_k)\right)_{1\leq t\leq n \atop 1\leq k\leq K}.
\end{align}
Note that under assumption C0, for all $k\in\wb{K}$,
\begin{align*}
	\big|(\widehat{\eta}_k^{(t)}-\widebar{\widehat{\eta}}_k)-(\eta_k^{(t)}-\mu_k)\big|
	\leq
	\big|\widehat{\eta}_k^{(t)}-\eta_k^{(t)}\big| + \big|\widebar{\widehat{\eta}}_k-\mu_k\big|
	\leq
	2\delta+\big|\widebar{\eta}_k-\mu_k\big|,
\end{align*}
which means
\begin{align} \label{ineq: norm_tildeR}
	|\widetilde{\bm{R}}|_\infty\leq 2\delta + \max_{k\in\wb{K}}|\widebar{\eta}_k-\mu_k|.
\end{align}
Together with \eqref{Pineq: maxMu}, we have
\begin{align} \label{Pineq: norm_tildeR}
	P\big(|\widetilde{\bm{R}}|_\infty\leq 2\delta + \delta_\mu\big)
	\geq
	1 - \sqrt{\frac{2}{\pi}}\frac{K\sigma}{\sqrt{n}\delta_\mu}\exp\left( -\frac{n\delta_\mu^2}{2\sigma^2}\right).
\end{align}
Using \eqref{Pineq: norm_tildeR}, Corollary~\ref{corollary: mb0_asy} (or Theorem~\ref{thm: ds0_asy}) and Lemma~\ref{lemma: net_etaDiff}, we have Corollary~\ref{corollary: net_mb_asy} and Corollary~\ref{corollary: net_ds_asy}.