# Distributed Sparse Linear Regression with Sublinear Communication

Chen Amiraz, Robert Krauthgamer and Boaz Nadler
Weizmann Institute of Science

###### Abstract

We study the problem of high-dimensional sparse linear regression in a
distributed setting under both computational and communication constraints.
Specifically, we consider a star topology network whereby several machines are
connected to a fusion center, with whom they can exchange relatively short
messages. Each machine holds noisy samples from a linear regression model with
the same unknown sparse $d$-dimensional vector of regression coefficients
$\bm{\theta}$. The goal of the fusion center is to estimate the vector
$\bm{\theta}$ and its support using few computations and limited communication
at each machine. In this work, we consider distributed algorithms based on
Orthogonal Matching Pursuit (OMP) and theoretically study their ability to
exactly recover the support of $\bm{\theta}$. We prove that under certain
conditions, even at low signal-to-noise-ratios where individual machines are
unable to detect the support of $\bm{\theta}$, distributed-OMP methods
correctly recover it with total communication sublinear in $d$. In addition,
we present simulations that illustrate the performance of distributed OMP-
based algorithms and show that they perform similarly to more sophisticated
and computationally intensive methods, and in some cases even outperform them.

##  1 Introduction

Sparse linear regression and feature selection are common problems in machine
learning, statistics and signal processing [21]. A common task is to identify,
out of a large set of explanatory variables, a small subset which is most
informative for a given prediction problem.

In this work we focus on the popular sparse linear regression model,

$y=\mathbf{x}^{T}\bm{\theta}+\sigma\xi,$ (1)

where $y\in\mathbb{R}$ is the response, $\mathbf{x}\in\mathbb{R}^{d}$ is a
vector of explanatory variables, $\bm{\theta}\in\mathbb{R}^{d}$ is the unknown
vector of regression coefficients, $\xi\in\mathbb{R}$ is a standard normal
random variable, i.e., $\xi\sim\mathcal{N}\left(0,1\right)$, and $\sigma>0$ is
the noise level. The coefficients vector $\bm{\theta}$ is assumed to be high-
dimensional, namely $d\gg 1$, and yet sparse, with sparsity level
$K=\left\lVert\bm{\theta}\right\rVert_{0}\ll d$. In a centralized setting,
given $N$ i.i.d. samples $\left(\mathbf{x}_{i},y_{i}\right)$ from (1), the
goal is to accurately estimate $\bm{\theta}$ as well as its support
$S=\\{i\mid\theta_{i}\neq 0\\}$. Over the past few decades many different
methods have been proposed and analyzed to solve this problem, including
combinatorial algorithms, greedy approaches and regularization schemes [33,
46, 12, 35, 11, 48, 10, 37, 16, 8, 22, 2].

With recent advances in data acquisition and storage, in various applications
the collected samples are stored across different machines. Moreover, due to
communication or privacy constraints, the machines cannot send their samples
to other machines in the network. This directly leads to various distributed
learning and inference problems, see for example [52, 23] and references
therein.

A common distributed setting for sparse linear regression, which we also
consider in this work, is that there are $M$ machines connected in a star
topology to a fusion center, with each machine having $n$ i.i.d. samples from
the model (1). Some distributed methods attempt to recover the centralized
solution, that would have been computed by the fusion center, if it had access
to all $nM$ samples of the $M$ machines. Examples include optimization-based
methods [34, 28, 29, 18, 36], Bayesian methods [32, 24], and greedy methods
[44, 27, 39, 20, 15]. These methods are iterative and in general communication
intensive, as they may require many rounds to converge. In contrast, several
authors proposed schemes with few rounds, possibly just one, and showed that
for a wide range of problems, the resulting estimator has a risk comparable to
that of the centralized solution [40, 49, 23, 30]. Specifically, in the
context of the sparse linear regression model (1), [26] and [6] proposed a
single-round debiased-Lasso approach. They further proved that under suitable
conditions and for various norms, the resulting estimator achieves the same
error rate as that of the centralized solution. Yet, these methods have two
limitations: (i) the communication per machine is at least linear in $d$; and
(ii) the required computational resources are non-negligible, since to compute
the debiased-Lasso estimator each machine has to solve $d+1$ separate Lasso
problems. Recently, [5] proposed a debiased-Lasso based distributed approach
with much lower communication. Their idea is for each machine to apply a
threshold to its debiased-Lasso estimate, and send to the fusion center only
the indices that exceed the threshold.

We focus on distributed estimation of the sparse vector $\bm{\theta}$ in the
model (1), under the following setting: The $M$ end machines have both limited
processing power and a restricted communication budget. This setting is
motivated by modern applications which involve end machines that are
relatively computationally weak, and yet are required to jointly learn from
high-dimensional data, namely $d\gg 1$. The first constraint implies that
computationally intensive methods such as debiased Lasso may be infeasible or
prohibitively slow. The second constraint implies that many of the previously
mentioned methods are not applicable, as their communication per machine is at
least linear in $d$. Since the quantity of interest $\bm{\theta}$ is
$K$-sparse with $K\ll d$, this setting raises the following challenge: develop
a scheme that accurately estimates the regression vector $\bm{\theta}$ with
number of operations per machine linear in $d$ and communication cost
_sublinear_ in $d$, and derive theoretical guarantees for it.

A natural base algorithm to address the above challenge is Orthogonal Matching
Pursuit (OMP). Indeed, [17, 50, 45] proposed to estimate the support of
$\bm{\theta}$ via various distributed-OMP schemes, which are both
computationally fast and require little communication. In terms of theory, to
the best of our knowledge, the only work to derive support recovery guarantees
for distributed-OMP methods is [51]. However, their analysis was restricted to
a noise-less compressed-sensing setting, where the design matrices are all
random and independent across machines. In this case, there is an underlying
symmetry between all non-support variables, which greatly simplifies their
proof. In contrast, we consider deterministic design matrices, for which their
proof technique is not applicable.

In this paper we make the following contributions. On the theoretical front,
we prove in Section 4 that at low signal-to-noise ratios, where individual
machines fail to recover the support of $\bm{\theta}$, a certain distributed-
OMP scheme can nonetheless correctly recover it with communication sublinear
in $d$. Moreover, its total communication cost decreases _exponentially_ as
the signal-to-noise ratio increases. Our theoretical result holds for
deterministic (low coherence) design matrices and noisy observations, as in
(1). The primary theoretical challenge in analyzing deterministic design
matrices is that the matrices may be similar (or even the same) across
machines. This implies that there may be high correlations between the same
support variables and non-support variables in each machine. At low signal-to-
noise ratios, this may lead several machines to send the _same_ incorrect
support variable to the fusion center. Deriving a theoretical guarantee in
this case requires a different and more delicate analysis than that of
previous works. Specifically, to bound the probability that a non-support
variable is sent to the fusion center we use recently derived lower bounds on
the maximum of correlated Gaussian random variables [31]. While our analysis
makes some simplifying assumptions, it provides insight into why distributed-
OMP methods can achieve exact support recovery even at low SNR where
individual machines fail to do so. To complement our theoretical results, in
Section 5 we present simulations that compare the support-recovery success of
distributed-OMP algorithms to several other methods, including OMP on a single
machine, debiased Lasso [26, 6], and a thresholded variant [5]. As expected,
debiased Lasso achieves exact support recovery at lower SNR, at the expense of
significantly higher communication and computational costs. Interestingly, in
comparison to a communication-restricted thresholded variant of debiased
Lasso, distributed-OMP methods perform comparably, and in some cases even
outperforms it, while being significantly faster. We conclude the paper with a
discussion in Section 6.

#### Notation

We use the standard $O(\cdot),\Omega(\cdot),\Theta(\cdot)$ notation to hide
constants independent of the problem parameters and the notation
$\tilde{O}(\cdot)$ to hide terms that are at most polylogarithmic in $d$. For
functions $f,g$ of the dimension $d$, the notations $f=o(g)$ and $f\ll g$ mean
that $f(d)/g(d)\to 0$ as $d\to\infty$. The term exact recovery of the support
$S$ with high probability means that an estimator $\hat{S}$ satisfies
$\Pr\left[\hat{S}=S\right]\to 1$ as both $d\to\infty$ and the number of
machines $M=M(d)\to\infty$ at a suitable rate, as detailed in each theorem.
The smallest integer larger than or equal to $x$ is denoted $\lceil x\rceil$.
The complement of the standard normal cumulative distribution function is
denoted by $\Phi^{c}(t)=\Pr\left[Z>t\right]$ where $Z\sim\mathcal{N}(0,1)$.

##  2 Problem setup

As mentioned in the introduction, we consider linear regression with a sparse
coefficient vector in a distributed setting, where $M$ machines are connected
in a star topology to a fusion center. Each machine $m$ holds $n$ samples of
the sparse linear regression model (1), i.e., a design matrix
$\mathbf{X}^{(m)}\in\mathbb{R}^{n\times d}$ and a response vector
$\mathbf{y}^{(m)}\in\mathbb{R}^{n}$. These two quantities are hence linearly
related according to

$\mathbf{y}^{(m)}=\mathbf{X}^{(m)}\bm{\theta}+\sigma\bm{\xi}^{(m)},$ (2)

where $\bm{\xi}^{(m)}\in\mathbb{R}^{n}$ is a random Gaussian noise vector
$\bm{\xi}^{(m)}\sim\mathcal{N}\left(0,\mathbf{I}_{n}\right)$. While different
machines may have the same or similar design matrices, the noise vectors in
different machines are assumed to be independent. We denote the support of
$\bm{\theta}$ by $S=\operatorname{supp}(\bm{\theta})$ and assume it is _$K$
-sparse_, namely
$\left\lVert\bm{\theta}\right\rVert_{0}=\left|\operatorname{supp}(\bm{\theta})\right|=K$.

In our work we view the $M$ matrices
$\left\\{\mathbf{X}^{(m)}\right\\}_{m\in[M]}$ as deterministic. In the model
above, each column of the design matrix corresponds to a different explanatory
variable, and thus can have a different second moment. We address this by
appropriately scaling the matrix columns in the algorithms. That being said, a
good example to keep in the mind is the case where all underlying variables
are scaled to have second moment equal to one, so that all $d$ columns of all
$M$ matrices have norm approximately equal to $\sqrt{n}$. This holds, for
example, under a generative model where the rows of $\mathbf{X}^{(m)}$ are
drawn from a Gaussian distribution $\mathcal{N}(0,\Sigma)$ with all diagonal
entries $\Sigma_{ii}=1$.

The problem we consider is exact recovery of the support of $\bm{\theta}$,
under the constraints that the $M$ machines have limited computational
resources and a limited communication bandwidth with the fusion center. For
simplicity, we assume that the sparsity level $K$ is known to the center.

##  3 Distributed OMP-based algorithms

We consider two distributed schemes to estimate the support of $\bm{\theta}$
that are based on Orthogonal Matching Pursuit (OMP), also known as Forward-
Stepwise Regression [13, 38, 33]. OMP is a natural candidate for machines with
low computational resources, since it is one of the simplest and fastest
greedy methods for sparse recovery.

The distributed OMP-based schemes outlined below use the following subroutine
which performs a single step of the OMP algorithm: Given a matrix
$\mathbf{X}$, a vector $\mathbf{y}$, and a current support estimate $\hat{S}$,
the subroutine estimates $\hat{\bm{\theta}}$ as its best approximation of
$\bm{\theta}$ on the support $\hat{S}$ and computes a corresponding residual
vector $\mathbf{r}$. Next, it calculates an index $j\in[d]$ whose column
$\mathbf{x}_{j}$ has maximal correlation with this residual $\mathbf{r}$ and
outputs it. This subroutine, which we call OMP_Step, is outlined in Algorithm
1. For future use, note that an important property of OMP_Step is the
orthogonality of the residual to the columns in the estimated support set.
This implies that the output index of OMP_Step satisfies $j\notin\hat{S}$.

The simplest distributed OMP-based method is for each machine to separately
run OMP for $K$ steps and send its locally-computed indices to the fusion
center. The center then forms a global estimator by choosing the $K$ indices
that received the highest number of votes. Since we are interested in low-SNR
regimes where individual machines may fail to recover the support, we propose
a variant where each machine runs OMP for a larger number of steps and thus
outputs a support of size $L>K$. This scheme, which we call Distributed OMP
(D-OMP), is outlined in Algorithm 2.

1

input : matrix $\mathbf{X}\in\mathbb{R}^{n\times d}$, vector
$\mathbf{y}\in\mathbb{R}^{n}$, current support estimate $\hat{S}$

output : support index $j$

2

3calculate
$\hat{\bm{\theta}}=\arg\min\left\\{\left\|\mathbf{y}-\mathbf{X}\mathbf{z}\right\|_{2}:\mathbf{z}\in\mathbb{R}^{d},\operatorname{supp}(\mathbf{z})=\hat{S}\right\\}$

4 compute residual $\mathbf{r}=\mathbf{y}-\mathbf{X}\hat{\bm{\theta}}$

5 output index
$j=\arg\max\left\\{\frac{\lvert\langle\mathbf{x}_{i},\mathbf{r}\rangle\rvert}{\left\lVert\mathbf{x}_{i}\right\rVert}:i\in[d]\right\\}$

6

Algorithm 1 OMP_Step (perform a single step of OMP)

1

2At each machine $m=1,\dots,M$  __

input : matrix $\mathbf{X}^{(m)}\in\mathbb{R}^{n\times d}$, vector
$\mathbf{y}^{(m)}\in\mathbb{R}^{n}$, parameter $L$

output : message $\hat{S}^{(m)}_{L}$ to center

3 initialize $\hat{S}^{(m)}_{0}=\emptyset$

4 for  _round $t=1,\dots,L$_ do

5 calculate index $j^{(m,t)}=$OMP_Step
$\left(\mathbf{X}^{(m)},\mathbf{y}^{(m)},\hat{S}^{(m)}_{t-1}\right)$

6 update estimated support set
$\hat{S}^{(m)}_{t}=\hat{S}^{(m)}_{t-1}\cup\left\\{j^{(m,t)}\right\\}$

7 end for

8 send $\hat{S}^{(m)}_{L}$ to the center

9

10At the fusion center  __

input : messages $\left\\{\hat{S}^{(m)}_{L}\right\\}_{m\in[M]}$, sparsity
level $K$

output : estimated support $\hat{S}$

11 for each index $j\in\left[d\right]$, calculate the number of votes it
received
$\mathbf{v}_{j}=\sum_{m\in[M]}\mathbb{1}\left\\{j\in\hat{S}^{(m)}_{L}\right\\}$

12 compute permutation $\pi$ that sorts the indices by descending number of
votes, $\mathbf{v}_{\pi(1)}\geq\dots\geq\mathbf{v}_{\pi(d)}$

13 return $K$ indices with most votes
$\hat{S}=\left\\{\pi(1),\dots,\pi(K)\right\\}$

14

15

Algorithm 2 Distributed OMP (D-OMP)

1

2initialize $\hat{S}_{0}=\emptyset$

3

4for  _round $t=1,\dots,K$_ do

5 At each machine $m=1,\dots,M$  __

input : matrix $\mathbf{X}^{(m)}\in\mathbb{R}^{n\times d}$, vector
$\mathbf{y}^{(m)}\in\mathbb{R}^{n}$, estimated support set $\hat{S}_{t-1}$

output : message $j^{(m,t)}=$OMP_Step
$\left(\mathbf{X}^{(m)},\mathbf{y}^{(m)},\hat{S}_{t-1}\right)$ sent to the
center

6

7

8 At the fusion center  __

input : messages $\left\\{j^{(m,t)}\right\\}_{m\in[M]}$, sparsity level $K$

output : estimated support $\hat{S}_{t}$

9 for each index $j\in\left[d\right]$, calculate the number of votes it
received
$\mathbf{v}^{(t)}_{j}=\sum_{m\in[M]}\mathbb{1}\left\\{j=j^{(m,t)}\right\\}$

10 find index with most votes
$j_{t}=\arg\max\\{\mathbf{v}^{(t)}_{j}:j\in[d]\\}$

11 add $j_{t}$ to the current estimated support
$\hat{S}_{t}=\hat{S}_{t-1}\cup\\{j_{t}\\}$

12 if  _$t=K$_ then

13 output $\hat{S}_{K}$

14 else

15 send $j_{t}$ to all machines

16 end if

17

18

19

20 end for

Algorithm 3 Distributed Joint OMP (DJ-OMP)

A second distributed OMP-based approach, which we call Distributed Joint OMP
(DJ-OMP), uses $K$ rounds of communication to construct the estimated support
set one index at a time. Starting with an empty support set, at each
communication round $t=1,\dots,K$, the center sends the current support
estimate $\hat{S}_{t-1}$ to the $M$ machines. Subsequently, each machine calls
OMP_Step and sends the resulting index $j^{(m,t)}$ to the center. At the end
of each round, the center chooses an index that received the highest number of
votes, denoted $j_{t}$, and adds it to the estimated support set
$\hat{S}_{t}$. After $K$ rounds, the center outputs the final support estimate
$\hat{S}_{K}$. Since each call to OMP_Step yields an index not in the current
support set $\hat{S}_{t-1}$, at each round $t$ of DJ-OMP no machine votes for
an index that is already in the estimated support set $\hat{S}_{t-1}$, and
hence the center chooses a new index $j_{t}\notin\hat{S}_{t-1}$. This scheme
is outlined in Algorithm 3.

For future use, let us analyze the number of operations in a single execution
of OMP_Step. First, given a current support set $\hat{S}$, the calculation of
$\hat{\bm{\theta}}$ via least squares involves multiplying a
$\left|\hat{S}\right|\times n$ matrix by its transpose, and then inverting the
resulting $\left|\hat{S}\right|\times\left|\hat{S}\right|$ matrix. Second,
calculating the index $j$ most correlated to the residual requires $d$ inner
products of vectors in $\mathbb{R}^{n}$. When $\left|\hat{S}\right|$ is
sufficiently small, say $o\left(d^{1/3}\right)$, the computational cost of
OMP_Step is dominated by the latter step whose cost is $O(nd)$.

We now briefly compare the two schemes DJ-OMP and D-OMP with $L=K$. First, in
terms of computational complexity, in both schemes each machine performs the
same number of total operations. Thus, assuming that
$K=o\left(d^{1/3}\right)$, their computational complexity per machine is
$O(ndK)$. Second, in terms of communication, in both schemes each machine
sends a total of $K$ indices, and thus their total communication is
$\tilde{O}(KM)$. The main difference between the schemes is that D-OMP
requires only a single round of communication, whereas DJ-OMP uses $K$ rounds.
Hence, DJ-OMP requires synchronization and is slower in comparison to D-OMP.

Several previous works proposed distributed-OMP methods as well. For example,
[50] considered the same D-OMP scheme as we do, with $L=K$. In addition, their
DC-OMP algorithm is similar to DJ-OMP in the following sense. In DC-OMP, at
each round each machine performs a single OMP step and sends an index whose
column has maximal correlation with the residual. However, the fusion is
performed differently. Their fusion rule is to check if there are indices that
received at least two votes. If so, they are all taken into the estimated
support set. Otherwise, a random index that received one vote is added to the
estimated support set. In comparison, our fusion rule is based on majority
voting, such that only one index is added to the support set at each round.

A distributed-OMP approach for a slightly different setting where each machine
has its own regression vector $\bm{\theta}^{(m)}$ was proposed in [45]. The
support sets of the $M$ vectors $\bm{\theta}^{(m)}$ are assumed to be similar,
and the $M$ machines are connected in a general topology without a fusion
center. The size of the common support set $K$ is also assumed to be known to
the machines. In their scheme, at each round $t=1,\dots,K$, each machine
exchanges the entire estimated support set with its neighbors.

Other distributed methods, based on Subspace Pursuit, have also been proposed
[44, 27]. These greedy methods operate differently from distributed-OMP. At
each communication round, the entire support set is estimated by each machine
and sent to the center, which computes a fused estimate and sends it back to
the machines. This continues until convergence. The number of rounds is not
a-priori known and depends on the network topology and other problem
parameters. Empirically, for the parameter regimes considered in these papers,
the number of rounds until convergence was of the order of the sparsity level.
Thus, due to the increase in the per round communication cost, these methods
require a total communication cost higher by at least a factor of $O(K)$
compared to distributed-OMP.

Despite the simplicity of these distributed schemes, to the best of our
knowledge, only limited theoretical results have been derived for them. The
most relevant result is [51], who proved a support recovery guarantee for DC-
OMP. However, [51] considered a restricted noise-less compressed-sensing
setting, where the entries of the design matrices are all random and i.i.d.
across machines. In contrast, we consider a more general setting, where the
design matrices are deterministic and potentially structured, and the samples
are noisy. In the next section, we derive a support recovery guarantee for a
variant of DJ-OMP for this setting. Then, in section 5 we compare distributed-
OMP schemes with various competitors via simulations.

##  4 Theoretical results

In this section we derive a support recovery guarantee for a variant of DJ-
OMP. Specifically, we prove that if the non-zero entries of $\bm{\theta}$ are
sufficiently large in absolute value, then with high probability a variant of
DJ-OMP will recover the support set $S$ using communication sublinear in $d$.
For simplicity, we first analyze the case $K=1$. We prove in Theorem 1 that in
this case, under suitable conditions, DJ-OMP achieves exact support recovery.

Next, we consider the general case $K\geq 1$. Deriving theoretical guarantees
for DJ-OMP in this case is challenging, since at each machine the same data is
used at different rounds of DJ-OMP. This creates a statistical dependence
between the indices chosen in different rounds which significantly complicated
the analysis. Hence, in Theorem 2 we instead derive recovery guarantees for a
variant of DJ-OMP, denoted DJ-OMP∗, that at each round uses fresh machines
with statistically independent data. From a practical perspective, as we show
in Section 5, empirically DJ-OMP and DJ-OMP∗ have similar performance. All
proofs appear in Appendix A.

Towards stating our result, we first recall a few definitions relevant for
recovery guarantees of OMP, and discuss the interesting SNR regime, where
individual machines fail to recover the support, but a distributed algorithm
may still succeed to do so.

#### Distributed coherence condition

Recall that the coherence of a matrix $\mathbf{A}$ with columns
$\mathbf{a}_{1},\dots,\mathbf{a}_{d}$ is defined as

$\mu\left(\mathbf{A}\right)=\max_{i\neq
j}\frac{\left|\left\langle\mathbf{a}_{i},\mathbf{a}_{j}\right\rangle\right|}{\left\lVert\mathbf{a}_{i}\right\rVert_{2}\left\lVert\mathbf{a}_{j}\right\rVert_{2}}.$
(3)

A matrix $\mathbf{A}$ satisfies the Mutual Incoherence Property (MIP) with
respect to a sparsity level $K$ if

$\mu(\mathbf{A})<\frac{1}{2K-1}.$ (4)

A fundamental result by Tropp [47] is that in an ideal noise-free setting, the
MIP condition (4) is sufficient for exact support recovery by OMP.

In our distributed setting, each machine $m$ has its own design matrix
$\mathbf{X}^{(m)}$, with coherence $\mu^{(m)}=\mu(\mathbf{X}^{(m)})$. We
denote the maximal coherence of all $M$ machines by

$\mu_{\max}=\mu_{\max}(\mathbf{X}^{(1)},\dots,\mathbf{X}^{(m)})=\max_{m\in[M]}\mu^{(m)}.$
(5)

We say that a set of matrices $\mathbf{X}^{(1)},\dots,\mathbf{X}^{(m)}$
satisfies the max-MIP condition with respect to a sparsity level $K$ if

$\mu_{\max}<\frac{1}{2K-1}.$ (6)

Notice that if the max-MIP condition (6) holds, then the MIP condition (4)
holds at each of the machines. This in turn implies that in the absence of
noise, OMP at each machine will correctly recover the support of
$\bm{\theta}$. The coherence plays a key role for OMP recovery also in the
presence of noise, as we discuss next.

#### SNR regime

Next, we present a formal definition of the SNR in our distributed setting. We
then focus on an interesting regime, in which the SNR values are sufficiently
high for OMP to recover the support of $\bm{\theta}$ in a centralized setting
where the center has access to all the samples from all machines, and yet too
low for OMP at a single machine to individually recover it. For a $K$-sparse
vector $\bm{\theta}\in\mathbb{R}^{d}$, a matrix
$\mathbf{A}\in\mathbb{R}^{n\times d}$ with coherence $\mu$ and columns with
unit norm, let

$\theta_{{\mbox{\tiny crit}}}(\mu,n,d,K,\sigma)=\frac{\sigma\sqrt{2\log
d}}{1-(2K-1)\mu}.$ (7)

Notice that $\theta_{{\mbox{\tiny crit}}}$ is well defined under the MIP
condition (4).

In our quest to derive exact support recovery guarantees, we shall consider
vectors $\bm{\theta}$ whose non-zero entries are not arbitrarily small, but
rather lower bounded by $\theta_{\min}$, namely $\min_{k\in
S}\left|\theta_{k}\right|\geq\theta_{\min}$. For a matrix $\mathbf{A}$ with
unit column norms, we define the SNR as
$r=\left(\frac{\theta_{\min}}{\theta_{{\mbox{\tiny
crit}}}(\mu,n,d,K,\sigma)}\right)^{2}$. As shown in previous works, near the
value $r=1$, OMP (at a single machine) exhibits a phase transition from
failure to success of support recovery. If the SNR is slightly higher, i.e.,
$r>\left(1+\sqrt{\frac{\log K}{\log d}}\right)^{2}$, then w.h.p. OMP exactly
recovers the support $S$. In contrast, if the SNR is slightly lower, i.e.,
$r<\left(1-\sqrt{\frac{\log K}{\log d}}-\mu\right)^{2}$, then there exists a
matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$ with coherence $\mu$ and a
$K$-sparse $\bm{\theta}\in\mathbb{R}^{d}$ for which, given the design matrix
$\mathbf{A}$ and a response vector
$\mathbf{y}=\mathbf{A}\bm{\theta}+\sigma\bm{\xi}$, OMP fails w.h.p. to recover
the support of $\bm{\theta}$. In addition, this occurs empirically for several
common choices of $\mathbf{A}$ and $\bm{\theta}$, see [3] and earlier results
in [7].

In our distributed setting the matrices $\mathbf{X}^{(m)}$ are assumed to be
deterministic and do not necessarily have unit column norms. However, notice
that the model (2) is equivalent to the following model

$\mathbf{y}^{(m)}=\tilde{\mathbf{X}}^{(m)}\tilde{\bm{\theta}}^{(m)}+\sigma\bm{\xi}^{(m)},$
(8)

where each column $\tilde{\mathbf{x}}_{j}^{(m)}$ of the matrix
$\tilde{\mathbf{X}}^{(m)}$ is scaled to have unit norm, i.e.,
$\tilde{\mathbf{x}}_{j}^{(m)}=\frac{\mathbf{x}_{j}^{(m)}}{\left\lVert\mathbf{x}_{j}^{(m)}\right\rVert}$,
and accordingly
$\tilde{\theta}_{j}^{(m)}=\left\lVert\mathbf{x}_{j}^{(m)}\right\rVert\theta_{j}$.
Clearly, the support of each $\tilde{\bm{\theta}}^{(m)}$ is identical to that
of $\bm{\theta}$. We assume the following lower bound $\tilde{\theta}_{\min}$
on the non-zero entries of the scaled vectors $\tilde{\bm{\theta}}^{(m)}$

$\min_{m\in[M]}\min_{k\in
S}\left|\tilde{\theta}_{k}^{(m)}\right|\geq\tilde{\theta}_{\min}.$ (9)

Given the above discussion, in our distributed setting we define the SNR
parameter $r$ as follows

$r=\left(\frac{\tilde{\theta}_{\min}}{\theta_{{\mbox{\tiny
crit}}}(\mu_{\max},n,d,K,\sigma)}\right)^{2}.$ (10)

If $r>1$ then $\tilde{\theta}_{\min}>\theta_{{\mbox{\tiny
crit}}}(\mu^{(m)},n,d,K,\sigma)$ for every machine $m\in[M]$, and hence OMP in
any single machine would recover the support of $\bm{\theta}$ with high
probability.

Next, consider a centralized setting where all $N=Mn$ samples are available to
the fusion center. In this setting, one forms a response vector
$\mathbf{y}\in\mathbb{R}^{N}$ and measurement matrix
$\mathbf{X}\in\mathbb{R}^{N\times d}$ by column-stacking
$\mathbf{y}^{(1)},\dots,\mathbf{y}^{(M)}$ and
$\mathbf{X}^{(1)},\dots,\mathbf{X}^{(M)}$ respectively. If we were to
normalize each column of $\mathbf{X}$, then the corresponding lower bound on
the non-zero entries of the coefficients would be $\min_{k\in
S}\left\\{\left\lVert\tilde{\mathbf{x}}_{k}\right\rVert\left|\theta_{k}\right|\right\\}\geq\sqrt{M}\tilde{\theta}_{\min}$.
Thus, in this centralized setting OMP may fail if
$\sqrt{M}\tilde{\theta}_{\min}<\theta_{{\mbox{\tiny
crit}}}(\mu(\mathbf{X}),N,d,K,\sigma)$. Hence, an SNR regime that is
interesting to study in the distributed setting is

$\frac{1}{M}\left(\frac{1-(2K-1)\mu_{\max}}{1-(2K-1)\mu(\mathbf{X})}\right)^{2}<r<1.$
(11)

As we show next, for a non-trivial subrange of the SNR values given in Eq.
(11), exact support recovery by distributed-OMP algorithms is possible with
total communication sublinear in $d$. To this end, we consider slightly
different variants of D-OMP and DJ-OMP, where we assume that, in addition to
$K$, also $\mu_{\max}$, $\tilde{\theta}_{\min}$ and $\sigma$, and thus the SNR
parameter $r$ of Eq. (10), are known to the fusion center. The center can use
its knowledge of $r$ to calculate how many machines suffice for exact support
recovery with high probability. We denote this number by $M_{c}$. If the total
number of machines $M$ is sufficiently large, in particular $M\geq M_{c}$,
then the center can contact only a subset of size $M_{c}$ of the machines,
hence reducing the total communication cost.

######  Remark 1.

Several previous works considered distributed settings where each machine $m$
has a different coefficients vector $\bm{\theta}^{(m)}$, but they are all
supported on the same set $S$, see for example [17, 28, 29, 51, 27]. Notice
that our analysis can also cover this case, by replacing the assumption (9) on
the non-zero entries of $\bm{\theta}$ with the following assumption

$\min_{m\in[M]}\min_{k\in
S}\left\\{\left\lVert\mathbf{x}_{k}^{(m)}\right\rVert\left|\theta_{k}^{(m)}\right|\right\\}\geq\tilde{\theta}_{\min}.$

###  4.1 Support recovery guarantee for sparsity $K=1$

For simplicity, we first state a support recovery guarantee for the case
$K=1$. We use the notation $M_{c}^{(1)}$ for the number of contacted machines
in this case. Towards specifying it, we define the following quantity $F$,
which will be used in the proof to lower bound the probability that the
support variable is sent to the center by one of the machines,

$F\left(d,\mu,r\right)=\left(1-\frac{\sqrt{1-\mu^{2}}}{\sqrt{\pi\log
d}}d^{-\frac{\mu^{2}}{1-\mu^{2}}}\right)\Phi^{c}\left(\frac{1-\sqrt{r}}{1-\mu}\sqrt{2\log
d}\right).$ (12)

We can now define $M_{c}^{(1)}\left(d,\mu,r\right)$ as follows

$M_{c}^{(1)}\left(d,\mu,r\right)=\left\lceil\frac{8\log
d}{F\left(d,\mu,r\right)}\right\rceil.$ (13)

The logarithmic factor arises from an application of Chernoff bounds to obtain
a support recovery with high probability. In our analysis, we assume that
$d\gg 1$ and that $\mu_{\max}$ is small. This implies that

$M_{c}^{(1)}\left(d,\mu_{\max},r\right)=\tilde{O}\left(d^{\left(\frac{1-\sqrt{r}}{1-\mu_{\max}}\right)^{2}}\right).$
(14)

Hence, as long as $1-\sqrt{r}<1-\mu_{\max}$, or equivalently
$r>\mu_{\max}^{2}$, the total communication is sublinear in $d$.

To guarantee support recovery, we also need to upper bound the probability
that a non-support variable is sent to the fusion center. For this, we use a
recent result on the lower-tail probability of the maximum of correlated
Gaussian random variables [31]. As we describe later, these will depend on a
parameter $\epsilon=\epsilon(\mu_{\max})$, required to satisfy the following
condition

$\frac{\sqrt{\mu_{\max}}}{1+\sqrt{\mu_{\max}}}<\epsilon<1.$ (15)

Since $\mu_{\max}$ is small, $\epsilon$ can be chosen to be small as well.

Finally, we define a few quantities towards describing the relevant values for
the SNR $r$ that guarantee support recovery by DJ-OMP. Let

$Q_{0}\left(d,\mu\right)=\frac{\log\left(44\sqrt{2}\right)-\log\left(1-\frac{\sqrt{1-\mu^{2}}}{\sqrt{\pi\log
d}}d^{-\frac{\mu^{2}}{1-\mu^{2}}}\right)}{\log d},$ (16)

Note that $Q_{0}=O\left(\log^{-1}d\right)$. In addition, let

$Q_{1}\left(d,\mu,\epsilon\right)=\frac{1-\left(1-\mu\right)\left(\left(1-\epsilon\right)\sqrt{1-\mu}-\sqrt{Q_{0}}\right)}{1-2\mu},$
(17)

and

$Q_{2}\left(d,\mu\right)=\frac{\sqrt{2+2\mu}\left(1+\left(1-\mu\right)\sqrt{Q_{0}}\right)}{1-\mu+\sqrt{2+2\mu}}.$
(18)

As we shall see below, a sufficient condition for successful support recovery
by DJ-OMP is that
$r>\min\left\\{Q_{1}\left(d,\mu_{\max},\epsilon\right),Q_{2}\left(d,\mu_{\max}\right)\right\\}^{2}$.
Notice that if $d\gg 1$ and both $\mu_{\max},\epsilon\ll 1$, then
$Q_{1}\approx\epsilon$, and $Q_{2}\approx\frac{\sqrt{2}}{1+\sqrt{2}}$. Thus,
$\min\left\\{Q_{1},Q_{2}\right\\}<1$, which implies that with high probability
DJ-OMP recovers the support of vectors $\bm{\theta}$ with SNRs at which
individual machines fail to do so.

######  Theorem 1.

Suppose that $\mu_{\max}<\frac{1}{2}$, that $\epsilon$ satisfies Eq. (15), and
that

$r>\min\left\\{Q_{1}\left(d,\mu_{\max},\epsilon\right),Q_{2}\left(d,\mu_{\max}\right)\right\\}^{2}.$
(19)

In addition, suppose that $M\geq M_{c}^{(1)}\left(d,\mu_{\max},r\right)$.
Then, for sufficiently large $d=d(\epsilon)$, with probability at least
$1-2/d,$ DJ-OMP with $M_{c}^{(1)}$ machines recovers the support of a
$1$-sparse vector $\bm{\theta}$ using a single round and
$O\left(M_{c}^{(1)}\log d\right)$ total communication bits.

A few remarks are in place. First, note that when $K=1$, D-OMP with $L=1$
reduces to the same algorithm as DJ-OMP, and thus this result holds for it as
well. Second, since $K=1$, the condition $\mu_{\max}<\frac{1}{2}$ is more
restrictive than the max-MIP condition (6), but only by a constant factor.
Third, when $\mu_{\max}\ll 1$ condition (19) roughly translates to
$r\gtrsim\epsilon^{2}$, and hence $r\gtrsim\mu_{\max}$. Thus, there is a non-
trivial range of relatively low SNR values for which DJ-OMP succeeds to
recover the support, even though the probability of any single machine to do
so it very low. Furthermore, since $\mu_{\max}<1$, the condition
$r>\mu_{\max}$ also implies that $r>\mu_{\max}^{2}$. As discussed above, the
latter condition implies that not only does DJ-OMP succeed but also does so
with total communication sublinear in $d$.

###  4.2 Support recovery guarantee for general sparsity $K$

We now describe the changes required for stating the support recovery
guarantee for a general $K\geq 1$. First, to bound the effect of the
projection onto the subspace orthogonal to $\hat{S}$, we define the following
auxiliary quantities,

$N_{a}=1-\frac{\left(K-1\right)\mu_{\max}^{2}}{1-\left(K-2\right)\mu_{\max}},$
(20)

and

$N_{b}=1-\mu_{\max}^{2}-\frac{\left(K-1\right)\mu_{\max}^{2}\left(1+\mu_{\max}\right)^{2}}{1-\left(K-2\right)\mu_{\max}}.$
(21)

As we show in Section A.5, these terms are both strictly positive under the
max-MIP condition (6). Condition (15) on $\epsilon\left(\mu_{\max}\right)$ is
replaced by111In fact, as shown in the proof, it suffice to require that
$\frac{\sqrt{\mu_{\max}+1-N_{a}}}{1+\sqrt{\mu_{\max}+1-N_{a}}}<\epsilon<1$,
which depends explicitly on $K$, but for simplicity we use the looser
condition (22).

$\displaystyle\frac{\sqrt{2\mu_{\max}}}{1+\sqrt{2\mu_{\max}}}<\epsilon<1.$
(22)

The probability $F$ in Eq. (12) is replaced by

$F\left(d,K,\mu,r\right)=\left(1-\sqrt{\frac{{N_{b}}}{{\pi\log
d}}}d^{1-1/N_{b}}\right)\Phi^{c}\left(\frac{1-\sqrt{r}}{\sqrt{N_{a}}(1-\mu)}\sqrt{2\log
d}\right),$ (23)

The quantity $M_{c}^{(1)}$ in Eq. (13) is replaced by

$M_{c}^{(1)}\left(d,K,\mu,r\right)=\left\lceil\frac{8\log
d}{F\left(d,K,\mu,r\right)}\right\rceil,$ (24)

the term $Q_{0}$ of Eq. (16) is replaced by

$Q_{0}\left(d,K,\mu\right)=\frac{\log\left(44\sqrt{2}K\right)-\log\left(1-\sqrt{\frac{{N_{b}}}{{\pi\log
d}}}d^{1-1/N_{b}}\right)}{\log d},$ (25)

the term $Q_{1}$ of Eq. (17) is replaced by

$Q_{1}\left(d,K,\mu,\epsilon\right)=\frac{1-\left(1-\mu\right)\sqrt{N_{a}}\left(\left(1-\epsilon\right)\sqrt{1-\mu}-\sqrt{Q_{0}}\right)}{1-2\mu
K\sqrt{N_{a}}\frac{1-\mu}{1-\left(2K-1\right)\mu}},$ (26)

and the term $Q_{2}$ of Eq. (18) is replaced by

$Q_{2}\left(d,K,\mu\right)=\frac{\sqrt{2+4\mu}\left(1+\sqrt{N_{a}}\left(1-\mu\right)\sqrt{Q_{0}}\right)}{\sqrt{N_{a}}\left(1-\mu\right)+\sqrt{2+4\mu}}.$
(27)

We can now state our support recovery guarantee for a general sparsity level
$K$. Denote by DJ-OMP∗ a variant of the DJ-OMP algorithm where at each round
the center contacts a fresh set of $M_{c}$ machines, and sends them the
current estimated support set. We prove in the following theorem that under a
slightly stricter condition than the max-MIP condition (6) and for
sufficiently high SNR, the DJ-OMP∗ algorithm recovers the support of
$\bm{\theta}$ with high probability.

######  Theorem 2.

Suppose that $(4K-1)\mu_{\max}-2K\mu_{\max}^{2}<1$, that $\epsilon$ satisfies
Eq. (22), and that

$r>\min\left\\{Q_{1}\left(d,K,\mu_{\max},\epsilon\right),Q_{2}\left(d,K,\mu_{\max}\right)\right\\}^{2}.$
(28)

In addition, suppose that $M\geq
K^{2}M_{c}^{(1)}\left(d,K,\mu_{\max},r\right)$. Then, for sufficiently large
$d=d(\epsilon)$, with probability at least $1-2K/d,$ DJ-OMP∗ with
$M_{c}=KM_{c}^{(1)}\left(d,K,\mu_{\max},r\right)$ fresh machines in each round
recovers the support of a $K$-sparse vector $\bm{\theta}$ using $K$ rounds and
$O\left(K^{3}M_{c}^{(1)}\log d\right)$ total communication bits.

Notice that for the recovery guarantee for general $K$ above, the number of
contacted machines in each round is larger by a factor of $K$ compared to the
case of a single support variable. The reason is that the lower bound that we
derive is on the number of votes received for all the support variables
combined. To lower bound the number of votes for a single support variable
using a pigeon hole argument, the number of machines must increase
accordingly. However, if $\bm{\theta}$ is sufficiently sparse, for example if
$K$ is polylogarithmic in $d$, then the remarks that followed Theorem 1 are
still in place. Specifically, in this case the total communication is
sublinear in $d$ as well.

We now briefly compare our results with some related works. Recently, [4]
studied a related problem of distributed sparse mean estimation, which can be
viewed as a special case of distributed sparse linear regression with an
orthogonal measurement matrix. For this case, they showed that for a wide
range of SNR values, exact support recovery is possible with total
communication sublinear in $d$. However, their proofs rely on the assumption
that the matrix is orthogonal, and cannot be easily generalized to incoherent
matrices. As mentioned above, [51] considered a compressed sensing setting
whereby the matrices are incoherent, but assumed to be random with entries
drawn i.i.d. from the same distribution. In both of these papers, the
probability for the fusion center to select a non-support index is uniform,
which greatly simplifies the analysis. We show that even without this
uniformity, distributed-OMP algorithms can achieve exact support recovery with
sublinear communication.

##  5 Simulation results

Figure 1 compares several algorithms with different computation and
communication costs, as outlined in Table 1. The first algorithm simply
performs OMP on a single machine, ignoring all other $M-1$ machines (blue).
The second algorithm is Deb-Lasso where each machine computes a
$d$-dimensional debiased-Lasso estimate of $\bm{\theta}$ and sends it to the
center. Subsequently, the center computes the average of the $M$ vectors and
returns its top $K$ indices (red) [26, 6]. The third algorithm is Deb-Lasso-K,
a variant of the method of [5], whereby each machine sends the top $K$ indices
of its debiased-Lasso estimate (yellow). The distributed-OMP algorithms that
we compare, described in Section 3, are D-OMP with $L=K$ rounds (solid
purple), D-OMP with $L=2K$ rounds (dotted purple), DJ-OMP (solid green) and
DJ-OMP∗ (dotted green). In addition, we present a lower bound that corresponds
to $\theta_{{\mbox{\tiny crit}}}(\mu(\mathbf{X}),N,d,K,\sigma)$ defined in Eq.
(7), below which even in the centralized setting exact support recovery by OMP
fails with high probability (dashed black).

Algorithm |  | Communication cost |  | Main computational cost
---|---|---|---|---
Single OMP |  | $O\left(K\log d\right)$ |  |  $O\left(ndK\right)$ operations at a single machine
Deb-Lasso |  | $\tilde{O}\left(Md\right)$ |  | solving $d+1$ Lasso optimization problems at each machine
Deb-Lasso-K |  | $\tilde{O}\left(MK\right)$ |
D-OMP |  | $\tilde{O}\left(MK\right)$ |  | $O\left(ndK\right)$ operations at each machine
DJ-OMP |  |
DJ-OMP∗ |  |
Table 1: Algorithms compared in simulations and their respective communication
and computational costs. The computational cost of OMP-based algorithms is
under the assumption that $\bm{\theta}$ is sufficiently sparse, say $K\ll
d^{1/3}$. We further assume that the parameter $L$ of D-OMP is $L=O(K)$.

(a) Setting 1: Entries of each measurement matrix are drawn i.i.d.
($\alpha=0$)

(b) Setting 2: Columns of each measurement matrix are slightly correlated
($\alpha=0.1$)

Figure 1: Comparison of success probability of various algorithms as a
function of $\theta_{\min}$.

We consider two different settings of parameters, as described below. Each
matrix $X^{(m)}$ is generated as follows. Each row is drawn independently from
$\mathcal{N}(\mathbf{0},\Sigma)$, where $\Sigma$ is a Toeplitz matrix with
$\Sigma_{ii}=1$ and $\Sigma_{ij}=\alpha^{\left|i-j\right|}$ for $i\neq j$ and
some $\alpha\in[0,1)$. In Setting 1, $\alpha=0$, i.e., each entry of each
matrix $X^{(m)}$ is an i.i.d. standard Gaussian random variable. In Setting 2,
$\alpha=0.1$, i.e., the columns of each matrix $X^{(m)}$ are slightly
correlated, as described above. The vector $\bm{\theta}$ is generated with
sparsity level $K=3$, such that
$\bm{\theta}=\theta_{\min}\cdot[1,-2,3,0,\dots,0]^{T}$. In both settings the
dimension is $d=2000$, the number of samples per machine is $n=1800$, the
number of machines is $M=70$, the noise level is $\sigma=1$. DJ-OMP∗ used $70$
fresh machines at each round, and thus effectively used data from $70K=210$
machines. The tuning parameter $\lambda$ of the debiased-Lasso methods, which
multiplies the $\ell_{1}$ term of each of the $d+1$ Lasso objectives, was set
to $\lambda=2\sigma\sqrt{\frac{\log d}{n}}$.

Figure 1 compares the empirical success probability of the aforementioned
algorithms as a function of $\theta_{\min}$ in the two settings outlined
above. Formally, for each algorithm $A$,

$p_{\text{ success
}}^{A}(\theta_{\min})=\frac{1}{J}\sum_{j=1}^{J}\mathbb{1}\left\\{\hat{S}_{j}^{A}(\theta_{\min})=S\right\\},$

where $\hat{S}_{j}^{A}(\theta_{\min})$ is the support set estimated by
algorithm $A$, for noise realization $j$ and lower bound $\theta_{\min}$ on
the non-zero coefficients of $\bm{\theta}$, and $J$ is the total number of
noise realizations, set to $J=500$.

The simulation results reveal several interesting behaviors. First, as
expected, in all settings distributed-OMP algorithms outperform OMP in a
single machine, which has very low computational and communication costs, and
underperform compared to Deb-Lasso, which has very high computational and
communication costs. Second, all the distributed-OMP algorithms succeed at
non-trivial SNR values, where OMP in a single machine fails with high
probability. Third, both DJ-OMP and DJ-OMP∗ perform similarly, so our
simplifying assumption that fresh machines are used in each round is not
needed in practice. Fourth, DJ-OMP and DJ-OMP∗ perform similarly to D-OMP with
$L=K$, and in Setting 2 they even slightly outperform it, which shows that for
some settings there is an advantage for selecting one support index at a time.
However, if one-shot communication is required, then at the expense of
increasing the communication budget (in this example by a factor of $2$),
D-OMP with $L=2K$ steps performs surprisingly well. Finally, in Setting 1
where each of the columns of the matrices $\mathbf{X}^{(m)}$ is drawn i.i.d.
from the unit sphere, distributed-OMP algorithms perform similarly to the much
more computationally intensive Deb-Lasso-K. Perhaps surprisingly, in the more
realistic Setting 2 where the columns are slightly correlated, distributed-OMP
methods even outperform it.

##  6 Discussion

In general, different algorithms for distributed sparse linear regression can
be compared based on three criteria: (i) the SNR above which they succeed with
high probability, (ii) communication costs, and (iii) computational
requirements. In Figure 2 we illustrate how distributed-OMP algorithms compare
with other methods proposed in the literature. Each circle represents an
algorithm. The horizontal axis represents the SNR, the vertical axis
represents the communication cost (on a logarithmic scale), and the size of
the circle represents the computational cost (also on a logarithmic scale).
The blue circle corresponds to OMP on a single machine, the red circle
corresponds to Deb-Lasso, the yellow circle corresponds to Deb-Lasso-K and the
purple circle corresponds to the distributed-OMP algorithms that we consider
in this paper.

As seen in the figure, the SNR and communication criteria exhibit a tradeoff.
When the SNR is sufficiently high for a single machine to recover the support
of $\bm{\theta}$ on its own, for example by OMP, then the total communication
cost of sending it to the center is only $K\log d$ bits. Notice that even in a
noise-less setting, for the fusion center to obtain the support, $K$ indices
must be sent to the center, so $K\log d$ bits is a lower bound on the required
communication. On the other hand, when the SNR is low, Deb-Lasso succeeds to
recover the support of $\bm{\theta}$ but incurs a communication cost of
$\tilde{O}\left(Md\right)$ bits, which is very intensive in a high-dimensional
setting. In general, at low-SNR values, we believe that no algorithm can
achieve exact support recovery with communication sublinear in $d$. Indeed,
for closely related problems, without a lower bound on the SNR, achieving the
centralized minimax $\ell_{2}$ risk or the centralized prediction error is not
possible with $o\left(d\right)$ communication [41, 43, 1]. However, for a wide
range of SNR values between these two extremes, we show that if the SNR is
known, then distributed-OMP algorithms do recover the support of $\bm{\theta}$
with total communication _sublinear_ in $d$. In fact, as implied by Eq. (14),
their total communication cost decreases _exponentially_ as the SNR increases.
An interesting open question it to determine the optimal rate at which the
communication cost decreases as a function of the SNR by any distributed
algorithm that achieves exact support recovery.

Computational constraints are important in settings where the machines have
limited processing power. As illustrated in Figure 2, distributed-OMP methods
and Deb-Lasso-K have comparable communication costs and they both succeed at
comparable SNR values. The major advantage of distributed-OMP methods over
Deb-Lasso-K is that they incur a much lower computational cost. In fact,
distributed-OMP algorithms can be represented by the dotted black curves that
connect the purple and blue circles in the figure – as the SNR increases, both
communication and computation costs decrease. The exact relationship between
communication costs and use of computational resources is another interesting
subject for future research.

Figure 2: Illustration of the tradeoff between SNR, communication, and
computation costs for several distributed sparse linear regression algorithms.
The distributed-OMP algorithms from Section 3 (in purple) strike a balance
between the SNRs at which they succeed, their communication and computational
costs.

## Acknowledgments

B.N. is the incumbent of the William Petschek professorial chair of
mathematics. R.K. is partially supported by ONR Award N00014-18-1-2364, by the
Israel Science Foundation grant #1086/18, and by a Minerva Foundation grant.
This research was partially supported by a grant from the Council for Higher
Education Competitive Program for Data Science Research Centers. We would like
to thank Rodney Fonseca for sharing with us his implementation of the Deb-
Lasso and Deb-Lasso-K algorithms.

## References

  * Acharya et al. [2019] Jayadev Acharya, Chris De Sa, Dylan Foster, and Karthik Sridharan.  Distributed learning with sublinear communication.  In _International Conference on Machine Learning_ , pages 40–50, 2019\.
  * Amir et al. [2021] Tal Amir, Ronen Basri, and Boaz Nadler.  The trimmed lasso: Sparse recovery guarantees and practical optimization by the generalized soft-min penalty.  _SIAM journal on mathematics of data science_ , 3(3):900–929, 2021.
  * Amiraz et al. [2020] Chen Amiraz, Robert Krauthgamer, and Boaz Nadler.  Tight recovery guarantees for orthogonal matching pursuit under gaussian noise.  _Information and Inference: A Journal of the IMA_ , 2020.
  * Amiraz et al. [2022] Chen Amiraz, Robert Krauthgamer, and Boaz Nadler.  Distributed sparse normal means estimation with sublinear communication.  _Information and Inference: A Journal of the IMA_ , 2022.
  * Barghi et al. [2021] Hanie Barghi, Amir Najafi, and Seyed Abolfazl Motahari.  Distributed sparse feature selection in communication-restricted networks.  _arXiv preprint arXiv:2111.02802_ , 2021.
  * Battey et al. [2018] Heather Battey, Jianqing Fan, Han Liu, Junwei Lu, and Ziwei Zhu.  Distributed testing and estimation under sparse high dimensional models.  _Annals of statistics_ , 46(3):1352, 2018.
  * Ben-Haim et al. [2010] Zvika Ben-Haim, Yonina C Eldar, and Michael Elad.  Coherence-based performance guarantees for estimating a sparse vector under random noise.  _IEEE Transactions on Signal Processing_ , 58(10):5030–5043, 2010.
  * Bertsimas et al. [2016] Dimitris Bertsimas, Angela King, and Rahul Mazumder.  Best subset selection via a modern optimization lens.  _The annals of statistics_ , 44(2):813–852, 2016\.
  * Birnbaum [1942] Zygmunt Wilhelm Birnbaum.  An inequality for Mill’s ratio.  _The Annals of Mathematical Statistics_ , 13(2):245–246, 1942.
  * Blumensath and Davies [2008] Thomas Blumensath and Mike E Davies.  Iterative thresholding for sparse approximations.  _Journal of Fourier analysis and Applications_ , 14(5):629–654, 2008.
  * Candes and Tao [2005] Emmanuel J Candes and Terence Tao.  Decoding by linear programming.  _IEEE Transactions on Information Theory_ , 51(12):4203–4215, 2005.
  * Chen et al. [2001] Scott Shaobing Chen, David L Donoho, and Michael A Saunders.  Atomic decomposition by basis pursuit.  _SIAM review_ , 43(1):129–159, 2001.
  * Chen et al. [1989] Sheng Chen, Stephen A Billings, and Wan Luo.  Orthogonal least squares methods and their application to non-linear system identification.  _International Journal of Control_ , 50(5):1873–1896, 1989.
  * Chernoff [1952] Herman Chernoff.  A measure of asymptotic efficiency for tests of a hypothesis based on the sum of observations.  _The Annals of Mathematical Statistics_ , 23(4):493–507, 1952.
  * Chouvardas et al. [2015] Symeon Chouvardas, Gerasimos Mileounis, Nicholas Kalouptsidis, and Sergios Theodoridis.  Greedy sparsity-promoting algorithms for distributed learning.  _IEEE Transactions on Signal Processing_ , 63(6):1419–1432, 2015.
  * Dai and Milenkovic [2009] Wei Dai and Olgica Milenkovic.  Subspace pursuit for compressive sensing signal reconstruction.  _IEEE Transactions on Information Theory_ , 55(5):2230–2249, 2009.
  * Duarte et al. [2005] Marco F Duarte, Shriram Sarvotham, Dror Baron, Michael B Wakin, and Richard G Baraniuk.  Distributed compressed sensing of jointly sparse signals.  In _Asilomar Conference on Signals, Systems and Computers, 2005._ , pages 1537–1541. IEEE, 2005.
  * Fosson et al. [2016] Sophie M Fosson, Javier Matamoros, Carles Antón-Haro, and Enrico Magli.  Distributed recovery of jointly sparse signals under communication constraints.  _IEEE Transactions on Signal Processing_ , 64(13):3470–3482, 2016.
  * Gordon [1941] Robert D Gordon.  Values of Mills’ ratio of area to bounding ordinate and of the normal probability integral for large values of the argument.  _The Annals of Mathematical Statistics_ , 12(3):364–366, 1941.
  * Han et al. [2015] Puxiao Han, Ruixin Niu, and Yonina C Eldar.  Modified distributed iterative hard thresholding.  In _2015 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 3766–3770. IEEE, 2015.
  * Hastie et al. [2015] Trevor Hastie, Robert Tibshirani, and Martin Wainwright.  _Statistical learning with sparsity: the lasso and generalizations_.  CRC Press, 2015.
  * Hastie et al. [2020] Trevor Hastie, Robert Tibshirani, and Ryan Tibshirani.  Best subset, forward stepwise or lasso? analysis and recommendations based on extensive comparisons.  _Statistical Science_ , 35(4):579–592, 2020.
  * Jordan et al. [2019] Michael I Jordan, Jason D Lee, and Yun Yang.  Communication-efficient distributed statistical inference.  _Journal of the American Statistical Association_ , 114(526):668–681, 2019.
  * Khanna and Murthy [2016] Saurabh Khanna and Chandra R Murthy.  Decentralized joint-sparse signal recovery: A sparse bayesian learning approach.  _IEEE Transactions on Signal and Information Processing over Networks_ , 3(1):29–45, 2016.
  * Komatu [1955] Yûsaku Komatu.  Elementary inequalities for Mills’ ratio.  _Rep. Statist. Appl. Res. Un. Jap. Sci. Engrs_ , 4:69–70, 1955.
  * Lee et al. [2017] Jason D Lee, Qiang Liu, Yuekai Sun, and Jonathan E Taylor.  Communication-efficient sparse regression.  _The Journal of Machine Learning Research_ , 18(1):115–144, 2017.
  * Li et al. [2015] Gang Li, Thakshila Wimalajeewa, and Pramod K Varshney.  Decentralized and collaborative subspace pursuit: A communication-efficient algorithm for joint sparsity pattern recovery with sensor networks.  _IEEE Transactions on Signal Processing_ , 64(3):556–566, 2015.
  * Ling and Tian [2011] Qing Ling and Zhi Tian.  Decentralized support detection of multiple measurement vectors with joint sparsity.  In _2011 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 2996–2999. IEEE, 2011.
  * Ling et al. [2012] Qing Ling, Zaiwen Wen, and Wotao Yin.  Decentralized jointly sparse optimization by reweighted $\ell_{q}$ minimization.  _IEEE Transactions on Signal Processing_ , 61(5):1165–1170, 2012.
  * Liu et al. [2022] Zhan Liu, Xiaoluo Zhao, and Yingli Pan.  Communication-efficient distributed estimation for high-dimensional large-scale linear regression.  _Metrika_ , pages 1–31, 2022.
  * Lopes and Yao [2022] Miles E Lopes and Junwen Yao.  A sharp lower-tail bound for gaussian maxima with application to bootstrap methods in high dimensions.  _Electronic Journal of Statistics_ , 16(1):58–83, 2022.
  * Makhzani and Valaee [2013] Alireza Makhzani and Shahrokh Valaee.  Distributed spectrum sensing in cognitive radios via graphical models.  In _2013 5th IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP)_ , pages 376–379. IEEE, 2013\.
  * Mallat and Zhang [1993] Stéphane G Mallat and Zhifeng Zhang.  Matching pursuits with time-frequency dictionaries.  _IEEE Transactions on signal processing_ , 41(12):3397–3415, 1993.
  * Mateos et al. [2010] Gonzalo Mateos, Juan Andrés Bazerque, and Georgios B Giannakis.  Distributed sparse linear regression.  _IEEE Transactions on Signal Processing_ , 58(10):5262–5276, 2010.
  * Miller [2002] Alan Miller.  _Subset selection in regression_.  CRC Press, 2002.
  * Mota et al. [2011] João FC Mota, João MF Xavier, Pedro MQ Aguiar, and Markus Puschel.  Distributed basis pursuit.  _IEEE Transactions on Signal Processing_ , 60(4):1942–1956, 2011.
  * Needell and Tropp [2009] Deanna Needell and Joel A Tropp.  CoSaMP: Iterative signal recovery from incomplete and inaccurate samples.  _Applied and Computational Harmonic Analysis_ , 26(3):301–321, 2009.
  * Pati et al. [1993] Yagyensh Chandra Pati, Ramin Rezaiifar, and PS Krishnaprasad.  Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition.  In _Conference Record of The Twenty-Seventh Asilomar Conference on Signals, Systems and Computers_ , pages 40–44. IEEE, 1993.
  * Patterson et al. [2014] Stacy Patterson, Yonina C Eldar, and Idit Keidar.  Distributed compressed sensing for static and time-varying networks.  _IEEE Transactions on Signal Processing_ , 62(19):4931–4946, 2014.
  * Rosenblatt and Nadler [2016] Jonathan D Rosenblatt and Boaz Nadler.  On the optimality of averaging in distributed statistical learning.  _Information and Inference: A Journal of the IMA_ , 5(4):379–404, 2016.
  * Shamir [2014] Ohad Shamir.  Fundamental limits of online and distributed algorithms for statistical learning and estimation.  In _Advances in Neural Information Processing Systems_ , pages 163–171, 2014.
  * Šidák [1967] Zbyněk Šidák.  Rectangular confidence regions for the means of multivariate normal distributions.  _Journal of the American Statistical Association_ , 62(318):626–633, 1967.
  * Steinhardt and Duchi [2015] Jacob Steinhardt and John Duchi.  Minimax rates for memory-bounded sparse linear regression.  In _Conference on Learning Theory_ , pages 1564–1587, 2015.
  * Sundman et al. [2012] Dennis Sundman, Saikat Chatterjee, and Mikael Skoglund.  A greedy pursuit algorithm for distributed compressed sensing.  In _2012 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 2729–2732. IEEE, 2012.
  * Sundman et al. [2014] Dennis Sundman, Saikat Chatterjee, and Mikael Skoglund.  Distributed greedy pursuit algorithms.  _Signal Processing_ , 105:298–315, 2014.
  * Tibshirani [1996] Robert Tibshirani.  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society. Series B (Methodological)_ , 58:267–288, 1996.
  * Tropp [2004] Joel A Tropp.  Greed is good: Algorithmic results for sparse approximation.  _IEEE Transactions on Information Theory_ , 50(10):2231–2242, 2004.
  * Tropp and Gilbert [2007] Joel A Tropp and Anna C Gilbert.  Signal recovery from random measurements via orthogonal matching pursuit.  _IEEE Transactions on Information Theory_ , 53(12):4655–4666, 2007.
  * Wang et al. [2017] Jialei Wang, Mladen Kolar, Nathan Srebro, and Tong Zhang.  Efficient distributed learning with sparsity.  In _International Conference on Machine Learning_ , pages 3636–3645, 2017.
  * Wimalajeewa and Varshney [2013] Thakshila Wimalajeewa and Pramod K Varshney.  Cooperative sparsity pattern recovery in distributed networks via distributed-OMP.  In _2013 IEEE International Conference on Acoustics, Speech and Signal Processing_ , pages 5288–5292. IEEE, 2013.
  * Wimalajeewa and Varshney [2014] Thakshila Wimalajeewa and Pramod K Varshney.  OMP based joint sparsity pattern recovery under communication constraints.  _IEEE Transactions on Signal Processing_ , 62(19):5059–5072, 2014.
  * Wimalajeewa and Varshney [2017] Thakshila Wimalajeewa and Pramod K Varshney.  Application of compressive sensing techniques in distributed sensor networks: A survey.  _arXiv preprint arXiv:1709.10401_ , 2017.

##  Appendix A Proofs

In this appendix we prove the theorems stated in Section 4. For ease of
presentation, in Section A.1 we prove Theorem 1 which addresses the case
$K=1$. We then extend this proof to the general case $K\geq 1$ and prove
Theorem 2 in Section A.2. We conclude by proving various auxiliary lemmas in
Sections A.3-A.6.

Towards proving both theorems, we begin with analyzing the communication cost
of DJ-OMP∗, and then turn to proving that it succeeds to recover the support
of $\bm{\theta}$ with high probability. To this end, we present a few
preliminaries, state useful lemmas and outline the proof.

#### Communication cost of DJ-OMP∗

Recall that in round $t$, each of the $M_{c}=KM_{c}^{(1)}$ participating
machines sends exactly one index of size $\log d$ bits to the fusion center.
Then, the center sends the current estimated support set $\hat{S}_{t}$ of size
$t\log d$ bits to the next set of $KM_{c}^{(1)}$ machines. Thus,
$(t+1)KM_{c}^{(1)}\log d$ bits are transmitted at round $t$. After $K$ rounds,
the total communication of DJ-OMP∗ is $O(K^{3}M_{c}^{(1)}\log d)$ bits, as
stated in Theorem 2.

#### Preliminaries

We consider a single round of DJ-OMP∗, and hence we omit the notation $t$ from
now on. Assume that $\hat{S}\subset S$ at the beginning of the round. As
described in line 3 of Algorithm 3, each machine $m$ calls the subroutine
OMP_Step. First, it estimates a sparse vector supported on $\hat{S}$, i.e.,

$\hat{\bm{\theta}}^{(m)}=\arg\min_{\mathbf{z}\in\mathbb{R}^{d}}\left\|\mathbf{y}^{(m)}-\mathbf{X}^{(m)}\mathbf{z}\right\|_{2}\text{
s.t. }\operatorname{supp}(\mathbf{z})=\hat{S}.$ (29)

Then, it uses $\hat{\bm{\theta}}^{(m)}$ to calculate the residual vector

$\mathbf{r}^{(m)}=\mathbf{y}^{(m)}-\mathbf{X}^{(m)}\hat{\bm{\theta}}^{(m)}.$
(30)

Finally, the machine sends to the fusion center the index

$j^{(m)}=\arg\max_{i\in[d]}\lvert\langle\tilde{\mathbf{x}}^{(m)}_{i},\mathbf{r}^{(m)}\rangle\rvert,$
(31)

where
$\tilde{\mathbf{x}}^{(m)}_{i}=\frac{\mathbf{x}^{(m)}_{i}}{\left\lVert\mathbf{x}^{(m)}_{i}\right\rVert}$
is the $i$-th column of $\mathbf{X}^{(m)}$ divided by its norm.

As described in lines 5-6 of Algorithm 3, given the messages sent by all
$M_{c}$ machines, the fusion center sums the votes received for each index
$j\in[d]\setminus\hat{S}$ and computes

$j_{\text{center}}=\arg\max_{j\in[d]\setminus\hat{S}}\mathbf{v}_{j}.$

Towards proving that with high probability $j_{\text{center}}\in
S\setminus\hat{S}$, for every machine $m$ we define an additional SNR
parameter $\rho^{(m)}$ that corresponds to the current round. Denote

$\tilde{\theta}_{\max}^{(m)}=\max_{k\in
S\setminus\hat{S}}\left\\{\left\lVert\mathbf{x}_{k}^{(m)}\right\rVert\left|\theta_{k}\right|\right\\}.$
(32)

Similarly to the definition of $r$ in Eq. (10), we define

$\rho^{(m)}=\left(\frac{\tilde{\theta}_{\max}^{(m)}}{\theta_{{\mbox{\tiny
crit}}}\left(\mu_{\max},n,d,K,\sigma\right)}\right)^{2},$ (33)

where $\theta_{{\mbox{\tiny crit}}}$ is defined in Eq. (7). As discussed in
Section 4, if $\rho^{(m)}>1$, then machine $m$ can recover the support index
that corresponds to $\tilde{\theta}_{\max}^{(m)}$ with high probability [3].
Therefore, in what follows, we assume that $\rho^{(m)}\leq 1$ in all machines
$m\in[M]$.

#### Proof outline and lemmas

Next, we state several useful lemmas and provide a high level overview of the
proof. Suppose that at the start of one round of DJ-OMP∗, a subset
$\tilde{\cal{M}}$ of the machines of size $\left|\tilde{\cal{M}}\right|=M_{c}$
are contacted by the fusion center and given an input set $\hat{S}\subset S$.
For the center to choose an index $j_{\text{center}}\in S\setminus\hat{S}$, it
suffices that there exists some support index $k\in S\setminus\hat{S}$ that
received more votes than any non-support index, namely,

$\mathbf{v}_{k}>\max_{j\notin S}\mathbf{v}_{j}.$ (34)

A sufficient condition for (34) to occur is that for some threshold $t_{c}$,
both

$\mathbf{v}_{k}>t_{c},$ (35)

and

$\max_{j\notin S}\mathbf{v}_{j}<t_{c}.$ (36)

We will use the threshold

$t_{c}=\frac{\sum_{m\in\cal{M}}F\left(d,K,\mu_{\max},\rho^{(m)}\right)}{M_{c}^{(1)}\left(d,K,\mu_{\max},r\right)F\left(d,K,\mu_{\max},r\right)}4\log
d,$ (37)

where $F$ and $M_{c}^{(1)}$ are defined in Eqs. (23) and (24) respectively,
and $\cal{M}$ is a subset of the machines contacted by the fusion center of
size $\left|\cal{M}\right|=M_{c}^{(1)}\left(d,K,\mu_{\max},r\right)$. As we
show below, for this threshold both (35) and (36) hold with high probability.
Note that $\rho^{(m)}$ and $t_{c}$ are not assumed to be known to the center
and are only used in the proof.

The following Lemma 1 lower bounds the threshold $t_{c}$, which will be useful
in our proofs. Its proof follows directly from the definition of $F$ in Eq.
(23) and appears in Section A.3.

######  Lemma 1.

The threshold $t_{c}$ defined in Eq. (37) satisfies

$t_{c}\geq 4\log d.$ (38)

The following Lemma 2 states that if the expected number of votes for an index
$k\in S\setminus\hat{S}$ is sufficiently high, then event (35) occurs with
high probability. The next Lemma 3 shows that if the expected number of votes
for each non-support index $j\notin S$ is sufficiently low, then event (36)
occurs with high probability. The proofs of these lemmas are based on simple
applications of Chernoff bounds and appear in Section A.3 as well.

######  Lemma 2.

Assume that $\hat{S}\subset S$. Let $\mathbf{v}_{k}$ be the number of votes
for the index $k$ and let $t_{c}$ be the threshold in Eq. (37). If
$\textup{E}\left[\mathbf{v}_{k}\right]\geq 2t_{c}$, then

$\Pr\left[\mathbf{v}_{k}\leq t_{c}\right]\leq\frac{1}{d}.$

######  Lemma 3.

Assume that $\hat{S}\subset S$. Let $\mathbf{v}_{j}$ be the number of votes
for an index $j$ and let $t_{c}$ be the threshold in Eq. (37). If
$\textup{E}\left[\mathbf{v}_{j}\right]\leq\frac{t_{c}}{5}$ for all $j\notin
S$, then

$\Pr\left[\max_{j\notin S}\mathbf{v}_{j}\geq t_{c}\right]\leq\frac{1}{d}.$

It remains to bound $\textup{E}\left[\mathbf{v}_{j}\right]$ for $j\in
S\setminus\hat{S}$ and for $j\notin S$. Towards this goal, denote by
$p^{(m)}_{j}$ the probability that machine $m$ sent index $j$, namely

$p^{(m)}_{j}=\Pr\left[j^{(m)}=j\right],$ (39)

where $j^{(m)}$ is defined in Eq. (31).

Since $\textup{E}\left[\mathbf{v}_{j}\right]=\sum_{m}p_{j}^{(m)}$, it suffices
to bound the probability $p_{j}^{(m)}$. For ease of presentation, we first
derive these bounds for the case $K=1$ in Section A.1, and then extend them
for the general case $K\geq 1$ in Section A.2.

###  A.1 Proof of Theorem 1

When $K=1$, the initial estimated support set is $\hat{S}=\emptyset$ and thus
it trivially holds that $\hat{S}\subset S$. Moreover, OMP_Step simplifies to
the following procedure. At each contacted machine $m$, the residual is simply
the response vector, i.e., $\mathbf{r}^{(m)}=\mathbf{y}^{(m)}$. Thus, machine
$m$ sends to the fusion center the index

$j^{(m)}=\arg\max_{i\in[d]}\lvert\langle\tilde{\mathbf{x}}^{(m)}_{i},\mathbf{y}^{(m)}\rangle\rvert.$
(40)

Another simplification in the case $K=1$ is that the support set only contains
one index, which we denote by $k$, i.e., $S=\left\\{k\right\\}$. To conclude
the proof of Theorem 1, we derive a lower bound on the probability
$p_{k}^{(m)}$ for the support index $k$ in the following Lemma 4 and an upper
bound on the probability $p_{j}^{(m)}$ for each non-support index $j\notin S$
in the following Lemma 5. Their proofs appear in Section A.4.

######  Lemma 4.

Assume that $\left\|\bm{\theta}\right\|_{0}=K=1$ and let
$S=\left\\{k\right\\}=\operatorname{supp}\left\\{\bm{\theta}\right\\}$. For
each machine $m$ and the support index $k\in S$,

$p_{k}^{(m)}\geq F\left(d,\mu_{\max},\rho^{(m)}\right),$ (41)

where $p_{k}^{(m)}$ and $F$ are defined in Eqs. (39) and (12) respectively.

######  Lemma 5.

Assume that $\left\|\bm{\theta}\right\|_{0}=K=1$ and let
$S=\operatorname{supp}\left\\{\bm{\theta}\right\\}$. Further assume that
$\rho^{(m)}$ of Eq. (33) satisfies $\rho^{(m)}\leq 1$ for each machine $m$. If
$\epsilon$ satisfies Condition (15), the SNR parameter $r$ satisfies Condition
(19), and the dimension $d=d(\epsilon)$ is sufficiently large, then for each
machine $m$ and each non-support index $j\notin S$,

$p_{j}^{(m)}\leq\frac{F\left(d,\mu_{\max},\rho^{(m)}\right)}{11},$ (42)

where $p_{j}^{(m)}$ and $F$ are defined in Eqs. (39) and (12) respectively.

We now formally prove Theorem 1 by combining the above lemmas.

###### Proof of Theorem 1.

Let ${\cal M}$ be the set of $M_{c}^{(1)}\left(d,\mu_{\max},r\right)$
contacted machines. We first analyze the probability that event (35) occurs.
By Lemma 4, for the support index $k\in S$, its expected number of votes is
$\textup{E}\left[\mathbf{v}_{k}\right]=\sum_{m\in{\cal
M}}p_{k}^{(m)}\geq\sum_{m\in{\cal M}}F\left(d,\mu_{\max},\rho^{(m)}\right)$.
By the definitions of $t_{c}$ in Eq. (37) and
$M_{c}^{(1)}\left(d,\mu_{\max},r\right)$ in Eq. (13),

$\textup{E}\left[\mathbf{v}_{k}\right]\geq\left\lceil\frac{8\log
d}{F\left(d,\mu_{\max},r\right)}\right\rceil\frac{F\left(d,\mu_{\max},r\right)}{4\log
d}t_{c}\geq 2t_{c}.$

By Lemma 2, the event (35) occurs with probability at least $1-1/d$.

Next, we analyze the probability that event (36) occurs. Fix a non-support
index $j\notin S$. Since $\rho^{(m)}\leq 1$, then by Lemma 5, its expected
number of votes is $\textup{E}\left[\mathbf{v}_{j}\right]=\sum_{m\in{\cal
M}}p_{j}^{(m)}\leq\frac{1}{11}\sum_{m\in{\cal
M}}F\left(d,\mu_{\max},\rho^{(m)}\right)$. By the definitions of $t_{c}$ in
Eq. (37) and $M_{c}^{(1)}$ in Eq. (13),

$\textup{E}\left[\mathbf{v}_{j}\right]\leq\frac{1}{11}\left\lceil\frac{8\log
d}{F\left(d,\mu_{\max},r\right)}\right\rceil\frac{F\left(d,\mu_{\max},r\right)}{4\log
d}t_{c}<\frac{t_{c}}{5}.$

The last inequality is justified as follows. Recall that $\lceil x\rceil\leq
x+1$ for all $x$. Thus,

$\left\lceil\frac{8\log
d}{F\left(d,\mu_{\max},r\right)}\right\rceil\frac{F\left(d,\mu_{\max},r\right)}{4\log
d}\leq 2+\frac{F\left(d,\mu_{\max},r\right)}{4\log d}.$

Note that $F$ in Eq. (12) denotes a probability function and thus
$F\left(d,\mu_{\max},r\right)\leq 1$. Hence, when $d\geq 8$, then $\log d>2$,
and the term $\frac{F\left(d,\mu_{\max},r\right)}{4\log d}\leq\frac{1}{8}$.

By Lemma 3, the event (36) occurs with probability at least $1-1/d$. A union
bound completes the proof of Theorem 1. ∎

###  A.2 Proof of Theorem 2

We now turn to proving that DJ-OMP∗ succeeds to recover the support of
$\bm{\theta}$ with general sparsity level $K$ with high probability. The proof
relies on the following lemma, which bounds the probability that the center
chooses an incorrect index at a single round of the algorithm.

######  Lemma 6.

Let $\hat{S}$ be the support set given as input to a single round of DJ-OMP∗
and denote by $j_{\text{center}}$ the index chosen by the center at the end of
this round. Then, assuming that the conditions of Theorem 2 hold and that
$\hat{S}\subset S$, the index $j_{\text{center}}$ also belongs to the support
set $S$ with high probability. Specifically,

$\Pr\left[j_{\text{center}}\notin S\mid\hat{S}\subset S\right]\leq 2d^{-1}.$
(43)

We now show how Theorem 2 follows directly from Lemma 6.

###### Proof of Theorem 2.

Assume Lemma 6 indeed holds. Recall that DJ-OMP∗ starts with
$\hat{S}_{0}=\emptyset$, adds exactly one new index to the estimated support
set at each round, and runs for exactly $K$ rounds. Hence, failure to recover
the support $S$ can only occur if in at least one of the rounds the center
added a non-support index to $\hat{S}_{K}$. We decompose the event
$\hat{S}_{K}\neq S$ by conditioning on the first round $t$ at which the center
chose an incorrect index $j_{t}\notin S$. By the law of total probability,

$\displaystyle\Pr\left[\hat{S}_{K}\neq S\right]$ $\displaystyle=$
$\displaystyle\sum_{t=1}^{K}\Pr\left[j_{t}\notin S\mid\hat{S}_{t-1}\subset
S\right]\Pr\left[\hat{S}_{t-1}\subset S\right]$ $\displaystyle\leq$
$\displaystyle\sum_{t=1}^{K}\Pr\left[j_{t}\notin S\mid\hat{S}_{t-1}\subset
S\right].$

By Lemma 6, each of the terms above is smaller than $2d^{-1}$. Hence,
$\Pr\left[\hat{S}_{K}\neq S\right]\leq 2Kd^{-1}$, which completes the proof of
Theorem 2. ∎

We now prove Lemma 6. Since $\hat{S}\subset S$, we need to bound the
probability $p_{j}^{(m)}$ of Eq. (39) for $j\in S\setminus\hat{S}$ and for
$j\notin S$. We shall do so using the following two lemmas. The first one,
Lemma 7, in fact lower bounds a slightly different quantity $q^{(m)}$, which
denotes the probability that the index sent by machine $m$ belongs to the
support $S\setminus\hat{S}$, i.e.,

$q^{(m)}=\Pr\left[j^{(m)}\in S\setminus\hat{S}\right].$ (44)

Lemma 8 upper bounds $p_{j}^{(m)}$ for each $j\notin S$. Their proofs appear
in Section A.5.

######  Lemma 7.

For each machine $m$,

$q^{(m)}\geq F\left(d,K,\mu_{\max},\rho^{(m)}\right),$ (45)

where $q^{(m)}$ and $F$ are defined in Eqs. (44) and (23) respectively.

######  Lemma 8.

Assume that $\rho^{(m)}\leq 1$ of Eq. (33) satisfies $\rho^{(m)}\leq 1$ for
each machine $m$. If $\epsilon$ satisfies Condition (22), the SNR parameter
$r$ satisfies Condition (28), and the dimension $d=d(\epsilon)$ is
sufficiently large, then for each machine $m$ and each non-support index
$j\notin S$,

$p_{j}^{(m)}\leq\frac{F\left(d,\mu_{\max},K,\rho^{(m)}\right)}{11K}.$ (46)

where $p_{j}^{(m)}$ and $F$ are defined in Eqs. (39) and (23) respectively.

We now formally prove Lemma 6 by combining the above lemmas.

###### Proof of Lemma 6.

Let $\tilde{\cal{M}}$ be the set of contacted machines in the current round,
of size
$\left|\tilde{\cal{M}}\right|=KM_{c}^{(1)}\left(d,K,\mu_{\max},r\right)$. We
first analyze the probability that event (35) occurs. Since $\hat{S}\subset
S$, the undetected support set is $S\setminus\hat{S}$. Let
$\mathbf{v}(S\setminus\hat{S})=\sum_{k\in S\setminus\hat{S}}\mathbf{v}_{k}$ be
the total number of votes received for all undetected support indices
combined. By Lemma 7, the expected number of votes is
$\textup{E}\left[\mathbf{v}(S\setminus\hat{S})\right]\geq\sum_{m\in\tilde{\cal{M}}}q^{(m)}\geq\sum_{m\in\tilde{\cal{M}}}F\left(d,K,\mu_{\max},\rho^{(m)}\right)$.
By a pigeon hole argument, there exists a support index $k\in
S\setminus\hat{S}$ that received votes from a subset ${\cal
M}_{k}\subseteq\tilde{\cal{M}}$ of machines, of size at least
$M_{c}^{(1)}\left(d,K,\mu,r\right)$. By the definitions of $M_{c}^{(1)}$ in
Eq. (24) and $t_{c}$ in Eq. (37), for this support index $k$,

$\textup{E}\left[\mathbf{v}_{k}\right]\geq\sum_{m\in{\cal
M}_{k}}F\left(d,K,\mu_{\max},\rho^{(m)}\right)\geq\left\lceil\frac{8\log
d}{F\left(d,K,\mu_{\max},r\right)}\right\rceil\frac{F\left(d,K,\mu_{\max},\rho^{(m)}\right)}{4\log
d}t_{c}\geq 2t_{c}.$

By Lemma 2, the event (35) occurs with probability at least $1-1/d$.

Similarly to the proof of Theorem 1, Lemmas 3 and 8 imply that the event (36)
also occurs with probability at least $1-1/d$. The only change in the proof is
that the set of contacted machines is of size $KM_{c}^{(1)}$, but the factor
$K$ gets canceled out by the $1/K$ factor in Lemma 8. A union bound completes
the proof of Lemma 6. ∎

###  A.3 Proofs of Lemmas 1, 2 and 3

We first prove Lemma 1 and then use it to prove Lemmas 2 and 3.

###### Proof of Lemma 1.

Recall the definition of the threshold $t_{c}$ in Eq. (37). By the definition
of $F$ in Eq. (23), it can be rewritten as follows

$t_{c}=\frac{\sum_{m\in\cal{M}}\Phi^{c}\left(\frac{\sqrt{2\log
d}}{\sqrt{N_{a}}(1-\mu_{\max})}\left(1-\sqrt{\rho^{(m)}}\right)\right)}{M_{c}^{(1)}\left(d,\mu_{\max},r\right)\Phi^{c}\left(\frac{\sqrt{2\log
d}}{\sqrt{N_{a}}(1-\mu_{\max})}\left(1-\sqrt{r}\right)\right)}4\log d.$ (47)

Next, recall that by Eq. (32),
$\tilde{\theta}_{\min}\leq\tilde{\theta}_{\max}^{(m)}$, and thus
$r\leq\rho^{(m)}$, for each $m\in\cal{M}$. The bound (38) follows from
noticing that $\Phi^{c}(t)$ is monotonically decreasing in $t$. ∎

Lemmas 2 and 3 are slightly more general versions of [4, Lemmas A.3 and A.4].
Since in our setting the probabilities $p_{k}^{(m)}$ depend not only on the
index $k$ but also on the machine $m$, and since the threshold $t_{c}$ is
defined in a more general manner, we provide their proofs below for
completeness.

We will use the following Chernoff bounds.

######  Lemma 9 ([14]).

Suppose $X_{1},\dots,X_{d}$ are independent Bernoulli random variables and let
$X$ denote their sum. Then, for any $\delta\geq 0$,

$\Pr\left[X\geq\left(1+\delta\right)\textup{E}\left[X\right]\right]\leq
e^{-\frac{\delta^{2}\textup{E}\left[X\right]}{2+\delta}},$ (48)

and for any $0\leq\delta\leq 1$,

$\Pr\left[X\leq\left(1-\delta\right)\textup{E}\left[X\right]\right]\leq
e^{-\frac{\delta^{2}\textup{E}\left[X\right]}{2}}.$ (49)

Next, we first introduce a few notations. Denote the indicator that machine
$m$ sends index $k$ by $I^{(m)}_{k}=\mathbb{1}\left\\{j^{(m)}=k\right\\}$. The
number of votes that $k$ receives is thus
$\mathbf{v}_{k}=\sum_{m\in\cal{M}}I^{(m)}_{k}$. Further denote
$E_{k}=\textup{E}[\mathbf{v}_{k}]$.

Recall that the noises $\left\\{\bm{\xi}^{(m)}\right\\}_{m\in[M]}$ are
independent. Since $\hat{S}\subset S$ is fixed and DJ-OMP∗ uses fresh machines
at each round, the indicators $\left\\{I^{(m)}_{k}\right\\}_{m\in\cal{M}}$ are
independent of each other and of $\hat{S}$. We now combine Lemmas 9 and 1 to
prove Lemmas 2 and 3.

###### Proof of Lemma 2.

Let $\delta_{k}=1-\frac{t_{c}}{E_{k}}$. By the discussion above, we may apply
Lemma 9 to the sum $\mathbf{v}_{k}$. The Chernoff bound (49) implies that

$\Pr\left[\mathbf{v}_{k}<t_{c}\right]=\Pr\left[\mathbf{v}_{k}<(1-\delta_{k})E_{k}\right]\leq\exp\left(-\delta_{k}^{2}E_{k}/2\right).$

The assumption $E_{k}\geq 2t_{c}$ and Lemma 1 imply that

$\Pr\left[\mathbf{v}_{k}<t_{c}\right]\leq\Pr\left[\mathbf{v}_{k}<E_{k}/2\right]\leq\exp\left(-E_{k}/8\right)\leq\exp\left(-t_{c}/4\right)\leq
1/d.$

∎

###### Proof of Lemma 3.

Fix $j\notin S$ and let $\delta_{j}=\frac{t_{c}}{E_{j}}-1.$ The assumption
$E_{j}\leq\frac{t_{c}}{5}$ implies that $\delta_{j}\geq 4$, and thus
$\delta_{j}/(2+\delta_{j})\geq 2/3$. By the discussion above, we may apply
Lemma 9 to each sum $\mathbf{v}_{j}$. The Chernoff bound (48) implies that

$\Pr\left[\mathbf{v}_{j}>t_{c}\right]=\Pr\left[\mathbf{v}_{j}>(1+\delta_{j})E_{j}\right]\leq\exp\left(-\frac{\delta_{j}}{2+\delta_{j}}\delta_{j}E_{j}\right)\leq\exp\left(-\frac{8t_{c}}{15}\right).$

By Lemma 1, the above probability is smaller than $d^{-2}$. We conclude, by
applying a union bound, that

$\Pr\left[\max_{j\not\in
S}\mathbf{v}_{j}>t_{c}\right]\leq(d-K)\Pr\left[\mathbf{v}_{j}>t_{c}\right]\leq
1/d.$

∎

###  A.4 Proofs of Lemmas 4 and 5

We begin with a few definitions and notations. For a vector $\mathbf{u}$ and a
set of indices $\cal{R}$, let
$\mathbf{u}_{\mid\cal{R}}\in\mathbb{R}^{|\cal{R}|}$ be the restriction of the
vector $\mathbf{u}$ to $\cal{R}$. Similarly, for a matrix
$\mathbf{A}\in\mathbb{R}^{n\times d}$, let
$\mathbf{A}_{\mid\cal{R}}\in\mathbb{R}^{n\times|\cal{R}|}$ be a restriction of
the columns of the matrix $\mathbf{A}$ to $\cal{R}$. Further denote by
$\mathbf{A}^{\dagger}$ the Moore-Penrose pseudo inverse of the matrix
$\mathbf{A}$, i.e.,
$\mathbf{A}^{\dagger}=\left(\mathbf{A}^{T}\mathbf{A}\right)^{-1}\mathbf{A}^{T}$
and notice that $\mathbf{A}^{\dagger}\mathbf{A}=\mathbf{I}$. Lastly, recall
that $\tilde{\mathbf{X}}^{(m)}\in\mathbb{R}^{n\times d}$ is the scaled design
matrix in machine $m$ and denote by
$\mathbf{P}^{(m)}_{\mathcal{R}}\in\mathbb{R}^{n\times n}$ an orthogonal
projection onto the span of $\tilde{\mathbf{X}}^{(m)}_{\mid\mathcal{R}}$,
i.e.,

$\mathbf{P}^{(m)}_{\mathcal{R}}=\tilde{\mathbf{X}}^{(m)}_{\mid\mathcal{R}}\left(\tilde{\mathbf{X}}^{(m)}_{\mid\mathcal{R}}\right)^{\dagger}.$
(50)

For simplicity of notation, in Sections A.4-A.6 we fix a machine $m$ and thus
omit the index $m$ from the proofs.

In our proofs we shall use the following well known Lemma 10 on Gaussian tail
bounds, the technical Lemma 11, whose proof appears in Section A.6, and Lemma
12, which bounds the lower-tail probability of the maximum of correlated
normal random variables [31]. Lemma 12 is a new result which is used in the
proofs of Lemmas 5 and 8 to bound the probability that a non-support index is
sent to the center.

######  Lemma 10 (Gaussian tail bounds [19]).

For $t>0$,

$\frac{t}{\sqrt{2\pi}(t^{2}+1)}e^{-t^{2}/2}\leq\Phi^{c}(t)\leq\frac{1}{\sqrt{2\pi}t}e^{-t^{2}/2}.$
(51)

######  Lemma 11.

If $a,b\geq 0$, then

$\Phi^{c}\left(a+b\right)<\sqrt{2}e^{-b^{2}/2}\Phi^{c}\left(a\right).$

######  Lemma 12 ([31]).

Let
$\left(Z_{1},\dots,Z_{d}\right)\sim\mathcal{N}\left(\mathbf{0},\mathbf{\Sigma}\right)$
where $\mathbf{\Sigma}_{ii}=1$ for all $i\in[d]$ and
$\mathbf{\Sigma}_{ij}\leq\eta<1$ for some fixed $\eta>0$ for all $i\neq
j\in[d]$. Fix $\delta\in(0,1)$. There is a constant $C>0$ depending only on
$(\eta,\delta)$ such that

$\Pr\left[\max_{i\in[d]}Z_{i}<\delta\sqrt{2(1-\eta)\log d}\right]\leq
Cd^{-\frac{(1-\eta)(1-\delta)^{2}}{\eta}}(\log
d)^{\frac{1-\eta(2-\delta)-\delta}{2\eta}}.$ (52)

For simplicity of notation, we omit the dependence on the machine index $m$.
In Lemmas 4 and 5 the sparsity level of $\bm{\theta}$ is $K=1$. Let $k\in S$
be the support index of $\bm{\theta}$. By Eq. (8), in this case

$\mathbf{y}=\tilde{\theta}_{k}\tilde{\mathbf{x}}_{k}+\sigma\bm{\xi}.$ (53)

Recall that by its definition in Eq. (32),
$\tilde{\theta}_{\max}=\left\lVert\mathbf{x}_{k}\right\rVert\left|\theta_{k}\right|$.
By Eq. (33) for $\rho$ and Eq. (7) for $\theta_{{\mbox{\tiny crit}}}$ with
$K=1$,

$\tilde{\theta}_{\max}=\frac{\sigma\sqrt{2\rho\log d}}{1-\mu_{\max}}.$ (54)

We now prove the lemmas.

###### Proof of Lemma 4.

Recall that $p_{k}$, defined in Eq. (39), is the probability that the support
index $k$ is selected by OMP_Step. This occurs if $k$ has the highest
correlation with the response vector than that of any non-support index
$j\notin S$, i.e., $p_{k}$ is the probability that

$\left|\left\langle\tilde{\mathbf{x}}_{k},\mathbf{y}\right\rangle\right|\geq\max_{i\notin
S}\left|\left\langle\tilde{\mathbf{x}}_{i},\mathbf{y}\right\rangle\right|.$
(55)

To prove the lemma we need to lower bound the probability of event (55). We
decompose the noise $\bm{\xi}$ in Eq. (53) as the sum of two components, the
first
$\bm{\xi}_{\parallel}=\mathbf{P}_{k}\bm{\xi}=\left\langle\tilde{\mathbf{x}}_{k},\bm{\xi}\right\rangle\tilde{\mathbf{x}}_{k}$
is parallel to $\tilde{\mathbf{x}}_{k}$, namely
$\left\langle\tilde{\mathbf{x}}_{k},\bm{\xi}_{\parallel}\right\rangle=\left\langle\tilde{\mathbf{x}}_{k},\bm{\xi}\right\rangle$,
and the second
$\bm{\xi}_{\perp}=\bm{\xi}-\bm{\xi}_{\parallel}=\left(\mathbf{I}-\mathbf{P}_{k}\right)\bm{\xi}$,
is orthogonal to $\tilde{\mathbf{x}}_{k}$, i.e.,
$\left\langle\tilde{\mathbf{x}}_{k},\bm{\xi}_{\perp}\right\rangle=0$.

Next, we use this decomposition to bound each of the terms of event (55). The
LHS of (55) can be bounded by

$\displaystyle\left|\left\langle\tilde{\mathbf{x}}_{k},\mathbf{y}\right\rangle\right|$
$\displaystyle\geq$
$\displaystyle\operatorname{sign}\left(\tilde{\theta}_{k}\right)\left\langle\tilde{\mathbf{x}}_{k},\mathbf{y}\right\rangle=\operatorname{sign}\left(\tilde{\theta}_{k}\right)\left(\tilde{\theta}_{k}\left\langle\tilde{\mathbf{x}}_{k},\tilde{\mathbf{x}}_{k}\right\rangle+\sigma\left\langle\tilde{\mathbf{x}}_{k},\bm{\xi}\right\rangle\right)$
(56) $\displaystyle=$
$\displaystyle\tilde{\theta}_{\max}+\sigma\operatorname{sign}\left(\tilde{\theta}_{k}\right)\left\langle\tilde{\mathbf{x}}_{k},\bm{\xi}_{\parallel}\right\rangle,$

where the second step follows from Eq. (8), and the third step follows from
the definitions of $\tilde{\theta}_{\max}$ and $\bm{\xi}_{\parallel}$ above.

The RHS of (55) can be bounded by

$\displaystyle\max_{i\notin
S}\left|\left\langle\tilde{\mathbf{x}}_{i},\mathbf{y}\right\rangle\right|$
$\displaystyle=$ $\displaystyle\max_{i\notin
S}\left|\tilde{\theta}_{k}\left\langle\tilde{\mathbf{x}}_{i},\tilde{\mathbf{x}}_{k}\right\rangle+\sigma\left\langle\tilde{\mathbf{x}}_{i},\bm{\xi}_{\parallel}+\bm{\xi}_{\perp}\right\rangle\right|$
(57) $\displaystyle\leq$
$\displaystyle\left(\tilde{\theta}_{\max}+\sigma\left|\left\langle\tilde{\mathbf{x}}_{k},\bm{\xi}\right\rangle\right|\right)\max_{i\notin
S}\left\\{\left|\left\langle\tilde{\mathbf{x}}_{i},\tilde{\mathbf{x}}_{k}\right\rangle\right|\right\\}+\sigma\max_{i\notin
S}\left|\left\langle\tilde{\mathbf{x}}_{i},\bm{\xi}_{\perp}\right\rangle\right|$
$\displaystyle\leq$
$\displaystyle\left(\tilde{\theta}_{\max}+\sigma\left|\left\langle\tilde{\mathbf{x}}_{k},\bm{\xi}_{\parallel}\right\rangle\right|\right)\mu_{\max}+\sigma\max_{i\notin
S}\left|\left\langle\tilde{\mathbf{x}}_{i},\bm{\xi}_{\perp}\right\rangle\right|.$

where the first step follows from Eq. (8) and the decomposition of the noise,
the second step follows from the triangle inequality and the definitions of
$\tilde{\theta}_{\max}$ and $\bm{\xi}_{\parallel}$, and the last step follows
from the definition of $\mu_{\max}$. Combining Eq. (56) with Eq. (57) implies
that a sufficient condition for event (55) to hold is that

$\max_{i\notin
S}\left|\left\langle\tilde{\mathbf{x}}_{i},\bm{\xi}_{\perp}\right\rangle\right|\leq\operatorname{sign}\left(\tilde{\theta}_{k}\right)\left\langle\tilde{\mathbf{x}}_{k},\bm{\xi}_{\parallel}\right\rangle-\mu_{\max}\left|\left\langle\tilde{\mathbf{x}}_{k},\bm{\xi}_{\parallel}\right\rangle\right|+\frac{\tilde{\theta}_{\max}\left(1-\mu_{\max}\right)}{\sigma}.$

By Eq. (54), the previous event may be written as

$\max_{i\notin
S}\left|\left\langle\tilde{\mathbf{x}}_{i},\bm{\xi}_{\perp}\right\rangle\right|\leq\operatorname{sign}\left(\tilde{\theta}_{k}\right)\left\langle\tilde{\mathbf{x}}_{k},\bm{\xi}_{\parallel}\right\rangle-\mu_{\max}\left|\left\langle\tilde{\mathbf{x}}_{k},\bm{\xi}_{\parallel}\right\rangle\right|+\sqrt{2\rho\log
d}.$ (58)

For a pair of independent random variables $A,B$ and a threshold
$T\in\mathbb{R}$, it holds that

$\Pr\left[A\leq B\right]\geq\Pr\left[A\leq T\wedge B\geq
T\right]=\Pr\left[A\leq T\right]\cdot\Pr\left[B\geq T\right].$ (59)

Let $T=\sqrt{2\log d}$, denote the LHS of (58) by $A$ and its RHS by $B$ and
note that $A$ and $B$ are independent. Hence, it suffices to bound the
probabilities of $A\leq T$ and $B\geq T$.

We begin with bounding $\Pr\left[A\leq T\right]$. Fix $i\notin S$. By
definition $\bm{\xi}_{\perp}=\left(\mathbf{I}-\mathbf{P}_{k}\right)\bm{\xi}$.
Note that
$\left\langle\tilde{\mathbf{x}}_{i},\bm{\xi}_{\perp}\right\rangle=\left\langle\tilde{\mathbf{x}}_{i},\left(\mathbf{I}-\mathbf{P}_{k}\right)\bm{\xi}\right\rangle=\left\langle\left(\mathbf{I}-\mathbf{P}_{k}\right)\tilde{\mathbf{x}}_{i},\bm{\xi}\right\rangle.$
Normalizing the inner product results in a standard normal random variable
$Z_{i}=\frac{\left\langle\tilde{\mathbf{x}}_{i},\bm{\xi}_{\perp}\right\rangle}{\left\|\left(\mathbf{I}-\mathbf{P}_{k}\right)\tilde{\mathbf{x}}_{i}\right\|_{2}}\sim\mathcal{N}(0,1).$
By the definition of $\mu_{\max}$,

$\left\|\left(\mathbf{I}-\mathbf{P}_{k}\right)\tilde{\mathbf{x}}_{i}\right\|^{2}=\tilde{\mathbf{x}}_{i}^{T}\left(\mathbf{I}-\mathbf{P}_{k}\right)\tilde{\mathbf{x}}_{i}=1-\left\langle\tilde{\mathbf{x}}_{i},\tilde{\mathbf{x}}_{k}\right\rangle^{2}\geq\gamma_{1},$

where $\gamma_{1}=1-\mu_{\max}^{2}$. Hence,

$\Pr\left[A\leq T\right]\geq\Pr\left[\max_{i\notin
S}\left|Z_{i}\right|\leq\frac{T}{\gamma_{1}}\right].$

Since $\left\\{Z_{i}\right\\}_{i\notin S}$ are jointly Gaussian, by [42, Thm.
1],

$\Pr\left[\max_{i\notin
S}\left|Z_{i}\right|\leq\frac{T}{\gamma_{1}}\right]\geq\prod_{i\notin
S}\Pr\left[\left|Z_{i}\right|\leq\frac{T}{\gamma_{1}}\right].$

Using the Gaussian tail bound (51) with $T=\sqrt{2\log d}$,

$\Pr\left[\left|Z_{i}\right|\leq\frac{T}{\gamma_{1}}\right]\geq
1-\frac{\gamma_{1}}{\sqrt{\pi\log d}}d^{-\frac{1}{\gamma_{1}^{2}}}.$

Since $\left\\{Z_{i}\right\\}_{i\notin S}$ are identically distributed and by
Bernoulli’s inequality, i.e., $(1-a)^{d}\geq 1-da$ for all $d\geq 1$ and
$a\in[0,1]$,

$\Pr\left[A\leq T\right]\geq 1-\frac{\gamma_{1}}{\sqrt{\pi\log
d}}d^{1-\frac{1}{\gamma_{1}^{2}}}.$ (60)

We now bound $\Pr\left[B\geq T\right]$, where $B$ is the RHS of (58). By the
definition of $\bm{\xi}_{\parallel}$, since $\tilde{\mathbf{x}}_{k}$ has unit
norm, the inner product
$\left\langle\tilde{\mathbf{x}}_{k},\bm{\xi}_{\parallel}\right\rangle=\left\langle\tilde{\mathbf{x}}_{k},\bm{\xi}\right\rangle\sim\mathcal{N}\left(0,1\right)$.

By the law of total probability with $Z\sim\mathcal{N}\left(0,1\right)$,

$\displaystyle\Pr\left[B\geq T\right]$ $\displaystyle=$
$\displaystyle\Pr\left[\operatorname{sign}\left(\tilde{\theta}_{k}\right)\left\langle\tilde{\mathbf{x}}_{k},\bm{\xi}_{\parallel}\right\rangle-\mu_{\max}\left|\left\langle\tilde{\mathbf{x}}_{k},\bm{\xi}_{\parallel}\right\rangle\right|\geq
T-\sqrt{2\rho\log d}\right]$ $\displaystyle\geq$
$\displaystyle\Pr\left[\gamma_{2}\left|Z\right|\geq T-\sqrt{2\rho\log
d}\mid\operatorname{sign}\left(Z\right)=\operatorname{sign}\left(\tilde{\theta}_{k}\right)\right]\cdot\Pr\left[\operatorname{sign}\left(Z\right)=\operatorname{sign}\left(\tilde{\theta}_{k}\right)\right],$

where $\gamma_{2}=1-\mu_{\max}$. Since $Z$ is symmetric,
$\Pr\left[\operatorname{sign}\left(Z\right)=\operatorname{sign}\left(\tilde{\theta}_{k}\right)\right]=\frac{1}{2}$.
In addition, the magnitude of $Z$ is independent on its sign. Thus,

$\Pr\left[B\geq T\right]\geq\frac{1}{2}\Pr\left[\gamma_{2}\left|Z\right|\geq
T-\sqrt{2\rho\log d}\right].$

Recalling that $T=\sqrt{2\log d}$ yields

$\Pr\left[B\geq
T\right]\geq\Phi^{c}\left(\frac{1-\sqrt{\rho}}{\gamma_{2}}\sqrt{2\log
d}\right).$ (61)

Inserting the bounds (60) with $\gamma_{1}=1-\mu_{\max}^{2}$ and (61) with
$\gamma_{2}=1-\mu_{\max}$ into (59) completes the proof of Lemma 4. ∎

###### Proof of Lemma 5.

Fix a non-support index $j\notin S$. Recall that $p_{j}$, defined in Eq. (39),
is the probability that index $j$ is selected by OMP_Step. This occurs if $j$
has the highest correlation with the response vector, i.e.,

$p_{j}=\Pr\left[\left|\left\langle\tilde{\mathbf{x}}_{j},\mathbf{y}\right\rangle\right|>\max_{i\neq
j}\left|\left\langle\tilde{\mathbf{x}}_{i},\mathbf{y}\right\rangle\right|\right].$
(62)

In particular, for the $j$-th index to be chosen, the correlation of the
$j$-th column with the response vector must exceed both that of the support
column $k$, as well as that of any other non-support column
$i\notin\left\\{k,j\right\\}$. Indeed, in what follows we separately upper
bound

$\Pr\left[\left|\left\langle\tilde{\mathbf{x}}_{j},\mathbf{y}\right\rangle\right|>\max_{i\notin\left\\{k,j\right\\}}\left|\left\langle\tilde{\mathbf{x}}_{i},\mathbf{y}\right\rangle\right|\right]$
(63)

and

$\Pr\left[\left|\left\langle\tilde{\mathbf{x}}_{j},\mathbf{y}\right\rangle\right|>\left|\left\langle\tilde{\mathbf{x}}_{k},\mathbf{y}\right\rangle\right|\right],$
(64)

and then use the following inequality to upper bound (62) by their minimum.
Specifically, denote
$A=\left|\left\langle\tilde{\mathbf{x}}_{j},\mathbf{y}\right\rangle\right|$,
$B=\max_{i\notin\left\\{k,j\right\\}}\left|\left\langle\tilde{\mathbf{x}}_{i},\mathbf{y}\right\rangle\right|$
and
$C=\left|\left\langle\tilde{\mathbf{x}}_{k},\mathbf{y}\right\rangle\right|$,
then

$\Pr\left[A>\max\left\\{B,C\right\\}\right]\leq\min\left\\{\Pr\left[A>B\right],\Pr\left[A>C\right]\right\\}.$
(65)

For later use in both bounds, by the triangle inequality, the random variable
$A$ can be upper bounded as follows

$\left|\left\langle\tilde{\mathbf{x}}_{j},\mathbf{y}\right\rangle\right|=\left|\tilde{\theta}_{k}\left\langle\tilde{\mathbf{x}}_{j},\tilde{\mathbf{x}}_{k}\right\rangle+\sigma\left\langle\tilde{\mathbf{x}}_{j},\bm{\xi}\right\rangle\right|\leq\tilde{\theta}_{\max}\mu_{\max}+\sigma\left|\left\langle\tilde{\mathbf{x}}_{j},\bm{\xi}\right\rangle\right|.$
(66)

We now begin with (63). For each non-support index $i\notin S$ such that
$i\neq j$,

$\left|\left\langle\tilde{\mathbf{x}}_{i},\mathbf{y}\right\rangle\right|\geq\left\langle\tilde{\mathbf{x}}_{i},\mathbf{y}\right\rangle=\tilde{\theta}_{k}\left\langle\tilde{\mathbf{x}}_{i},\tilde{\mathbf{x}}_{k}\right\rangle+\sigma\left\langle\tilde{\mathbf{x}}_{i},\bm{\xi}\right\rangle\geq-\tilde{\theta}_{\max}\mu_{\max}+\sigma\left\langle\tilde{\mathbf{x}}_{i},\bm{\xi}\right\rangle.$

Combining the above bound with Eq. (66), rearranging the terms, and recalling
the relation between $\theta_{\max}$ and $\rho$ in (54) yields

$\Pr\left[\left|\left\langle\tilde{\mathbf{x}}_{j},\mathbf{y}\right\rangle\right|>\max_{i\notin\left\\{k,j\right\\}}\left|\left\langle\tilde{\mathbf{x}}_{i},\mathbf{y}\right\rangle\right|\right]\leq\Pr\left[\left|\left\langle\tilde{\mathbf{x}}_{j},\bm{\xi}\right\rangle\right|+2\mu_{\max}\frac{\sqrt{2\rho\log
d}}{1-\mu_{\max}}>\max_{i\notin\left\\{k,j\right\\}}\left\langle\tilde{\mathbf{x}}_{i},\bm{\xi}\right\rangle\right].$
(67)

Next, we use the following inequality which holds for any pair of random
variables $D,E$ and constant $T\in\mathbb{R}$,

$\Pr\left[D>E\right]\leq\Pr\left[D\geq T\right]+\Pr\left[E<T\right].$ (68)

By applying the above bound with
$T=(1-\epsilon)\sqrt{2\left(1-\mu_{\max}\right)\log d}$ for
$\epsilon\in\left(0,1\right)$ as in Eq. (15), we can upper bound (67) by

$\Pr\left[\left|\left\langle\tilde{\mathbf{x}}_{j},\bm{\xi}\right\rangle\right|\geq
a\sqrt{2\log
d}\right]+\Pr\left[\max_{i\notin\left\\{k,j\right\\}}\left\langle\tilde{\mathbf{x}}_{i},\bm{\xi}\right\rangle<(1-\epsilon)\sqrt{2\left(1-\mu_{\max}\right)\log
d}\right],$

where

$a=(1-\epsilon)\sqrt{1-\mu_{\max}}-\frac{2\mu_{\max}\sqrt{\rho}}{1-\mu_{\max}}.$

Since $\tilde{\mathbf{x}}_{j}$ has unit norm,
$\left\langle\tilde{\mathbf{x}}_{j},\bm{\xi}\right\rangle\sim\mathcal{N}\left(0,1\right)$.
Hence, the first term is bounded by

$2\Phi^{c}\left(a\sqrt{2\log d}\right).$ (69)

We now bound the second term. For the random variables
$\left\langle\tilde{\mathbf{x}}_{i},\bm{\xi}\right\rangle$, the $(i,l)$-th
entry of their covariance matrix is
$\textup{Cov}\left(\left\langle\tilde{\mathbf{x}}_{i},\bm{\xi}\right\rangle,\left\langle\tilde{\mathbf{x}}_{l},\bm{\xi}\right\rangle\right)=\left\langle\tilde{\mathbf{x}}_{i},\tilde{\mathbf{x}}_{l}\right\rangle$.
Thus, each of its diagonal entries is $1$ and each of its off-diagonal entries
is bounded in absolute value by $\mu_{\max}$. Hence, we can apply Lemma 12
with $\eta=\mu_{\max}$ and $\delta=1-\epsilon$, which gives the following
bound

$C\left(d-2\right)^{-\frac{1-\mu_{\max}}{\mu_{\max}}\epsilon^{2}}(\log(d-2))^{\frac{\epsilon-\mu_{\max}(1+\epsilon)}{2\mu_{\max}}}.$
(70)

We now show that (69) is larger than (70), and thus

$\Pr\left[\left|\left\langle\tilde{\mathbf{x}}_{j},\mathbf{y}\right\rangle\right|>\max_{i\notin\left\\{k,j\right\\}}\left|\left\langle\tilde{\mathbf{x}}_{i},\mathbf{y}\right\rangle\right|\right]\leq
4\Phi^{c}\left(a\sqrt{2\log d}\right).$ (71)

First note that if $\rho$ is sufficiently large such that $a\leq 0$, then (69)
is larger then 1, and thus larger than (70). Otherwise, $a>0$ and using the
lower bound for the Gaussian tail of (51), we may lower bound (69) by
$d^{-a^{2}-o(1)}$, where $o(1)$ hides factors that are asymptotically smaller
than 1. The term (70) can be upper bounded by $d^{-b^{2}+o(1)}$, where
$b=\sqrt{\frac{1-\mu_{\max}}{\mu_{\max}}}\epsilon$. Next, let us show that for
a fixed $\epsilon>0$, $b-a$ is bounded away from 0. This, in turn, implies
that for sufficiently large $d=d\left(\epsilon\right)$, (69) is larger than
(70). Indeed, under condition (15),
$\epsilon=\frac{\sqrt{\mu_{\max}}}{1+\sqrt{\mu_{\max}}}+\epsilon_{0}$ for some
$\epsilon_{0}>0$. Thus,
$b-a=\epsilon_{0}\sqrt{1-\mu_{\max}}\left(1+\frac{1}{\sqrt{\mu_{\max}}}\right)+\frac{2\mu_{\max}\sqrt{\rho}}{1-\mu_{\max}}$,
which is a sum of positive terms and hence bounded away from 0 as desired.
Therefore, condition (15) implies that (63) can be bounded by (71).

We now turn to bounding (64). For the support index $k$, by (53),

$\left|\left\langle\tilde{\mathbf{x}}_{k},\mathbf{y}\right\rangle\right|\geq\operatorname{sign}\left(\tilde{\theta}_{k}\right)\left\langle\tilde{\mathbf{x}}_{k},\mathbf{y}\right\rangle=\operatorname{sign}\left(\tilde{\theta}_{k}\right)\left(\tilde{\theta}_{k}\left\langle\tilde{\mathbf{x}}_{k},\tilde{\mathbf{x}}_{k}\right\rangle+\sigma\left\langle\tilde{\mathbf{x}}_{k},\bm{\xi}\right\rangle\right)=\tilde{\theta}_{\max}+\sigma\operatorname{sign}\left(\tilde{\theta}_{k}\right)\left\langle\tilde{\mathbf{x}}_{k},\bm{\xi}\right\rangle.$

Combining the bound above with (66) and inserting $\tilde{\theta}_{\max}$ in
Eq. (54), the probability (64) is upper bounded by

$\Pr\left[\left|\left\langle\tilde{\mathbf{x}}_{j},\bm{\xi}\right\rangle\right|-\operatorname{sign}\left(\tilde{\theta}_{k}\right)\left\langle\tilde{\mathbf{x}}_{k},\bm{\xi}\right\rangle>\sqrt{2\rho\log
d}\right].$ (72)

We now upper bound this probability. Let
$H=\left\langle\tilde{\mathbf{x}}_{j},\bm{\xi}\right\rangle$,
$G=\operatorname{sign}\left(\tilde{\theta}_{k}\right)\left\langle\tilde{\mathbf{x}}_{k},\bm{\xi}\right\rangle$
and $c=\sqrt{2\rho\log d}$.

For any pair of random variables $G,H$ and constant $c$,

$\Pr\left[\left|H\right|-G>c\right]\leq\Pr\left[H-G>c\right]+\Pr\left[-H-G>c\right].$
(73)

By their definition, $H,G$ are jointly Gaussian with mean zero and covariance
matrix

$\left(\begin{array}[]{cc}\sigma_{H}^{2}&\sigma_{HG}\\\
\sigma_{HG}&\sigma_{G}^{2}\end{array}\right).$

Hence, $H-G\sim\mathcal{N}(0,\sigma_{H}^{2}+\sigma_{G}^{2}-2\sigma_{HG})$ and
$-H-G\sim\mathcal{N}(0,\sigma_{H}^{2}+\sigma_{G}^{2}+2\sigma_{HG})$. Similarly
to the above discussion, the diagonal entries
$\sigma_{H}^{2}=\sigma_{G}^{2}=1$ and the off-diagonal entry
$\sigma_{HG}=\operatorname{sign}\left(\tilde{\theta}_{k}\right)\left\langle\tilde{\mathbf{x}}_{k},\tilde{\mathbf{x}}_{j}\right\rangle$.
Since $\left|\sigma_{HG}\right|\leq\mu_{\max}$, then by (73),

$\Pr\left[\left|H\right|-G>c\right]\leq\Phi^{c}\left(\frac{c}{\sqrt{2-2\sigma_{HG}}}\right)+\Phi^{c}\left(\frac{c}{\sqrt{2+2\sigma_{HG}}}\right)\leq
2\Phi^{c}\left(\frac{c}{\sqrt{2+2\mu_{\max}}}\right).$

Inserting $c=\sqrt{2\rho\log d}$ yields

$\Pr\left[\left|\left\langle\tilde{\mathbf{x}}_{j},\mathbf{y}\right\rangle\right|>\left|\left\langle\tilde{\mathbf{x}}_{k},\mathbf{y}\right\rangle\right|\right]\leq
2\Phi^{c}\left(\sqrt{\frac{2\rho\log d}{2+2\mu_{\max}}}\right).$ (74)

By Eq. (65), the probability (62) is at most the minimum between (71) and
(74). By the monotonicity of the Gaussian CDF, it is upper bounded by

$4\Phi^{c}\left(\max\left\\{\left((1-\epsilon)\sqrt{1-\mu_{\max}}-\frac{2\mu_{\max}\sqrt{\rho}}{1-\mu_{\max}}\right),\sqrt{\frac{\rho}{2+2\mu_{\max}}}\right\\}\sqrt{2\log
d}\right).$ (75)

Finally, to prove (42) of the lemma, we note that with $Q_{1}$ and $Q_{2}$
defined in Eqs. (17) and (18) respectively, condition (19) may be equivalently
written as

$\frac{1-\sqrt{r}}{1-\mu_{\max}}+\sqrt{Q_{0}}<\max\left\\{\left((1-\epsilon)\sqrt{1-\mu_{\max}}-\frac{2\mu_{\max}\sqrt{r}}{1-\mu_{\max}}\right),\sqrt{\frac{r}{2+2\mu_{\max}}}\right\\}.$
(76)

The definitions of $r$ and $\rho$ in Eqs. (10) and (33) imply that $\rho\geq
r$. Thus, $\rho$ satisfies condition (19) and hence condition (76). The RHS of
(76) is the same as the maximum in (75) above. Thus, (75) is upper bounded by

$4\Phi^{c}\left(\left(\frac{1-\sqrt{\rho}}{1-\mu_{\max}}+\sqrt{Q_{0}}\right)\sqrt{2\log
d}\right).$ (77)

Since $\rho\leq 1$, we can apply Lemma 11. Hence, by the definition of $Q_{0}$
in Eq. (16), and by the definition of $F$ in Eq. (12),

$\displaystyle p_{j}$ $\displaystyle\leq$ $\displaystyle
4\Phi^{c}\left(\left(\frac{1-\sqrt{\rho}}{1-\mu_{\max}}+\sqrt{Q_{0}}\right)\sqrt{2\log
d}\right)\leq
4\sqrt{2}d^{-Q_{0}}\Phi^{c}\left(\frac{1-\sqrt{\rho}}{1-\mu_{\max}}\sqrt{2\log
d}\right)$
$\displaystyle=4\sqrt{2}\frac{1-\frac{\sqrt{1-\mu^{2}}}{\sqrt{\pi\log
d}}d^{-\frac{\mu^{2}}{1-\mu^{2}}}}{44\sqrt{2}}\Phi^{c}\left(\frac{1-\sqrt{\rho}}{1-\mu_{\max}}\sqrt{2\log
d}\right)=\frac{F\left(d,\mu_{\max},\rho\right)}{11},$

which completes the proof of Lemma 5. ∎

###  A.5 Proof of Lemmas 7 and 8

We first make a few definitions and present a useful technical lemma. We begin
by rewriting the residual $\mathbf{r}^{(m)}$ using the notations introduced in
Section A.4. Recall that given an input support set $\hat{S}$, each machine
$m$ estimates its vector $\hat{\bm{\theta}}^{(m)}$ by solving the least
squares problem (29). Thus,
$\operatorname{supp}\left(\hat{\bm{\theta}}^{(m)}\right)=\hat{S}$ and

$\hat{\bm{\theta}}^{(m)}_{\mid\hat{S}}=\left(\mathbf{X}^{(m)}_{\mid\hat{S}}\right)^{\dagger}\mathbf{y}^{(m)}.$

Denote by $\tilde{\bm{\xi}}^{(m)}$ the projection of the noise
$\bm{\xi}^{(m)}$ to the subspace orthogonal to the span of the columns of
$\mathbf{X}^{(m)}_{\mid\hat{S}}$, i.e.,
$\tilde{\bm{\xi}}^{(m)}=\left(\mathbf{I}-\mathbf{P}_{\hat{S}}^{(m)}\right)\bm{\xi}^{(m)}$.
Given that $\hat{S}\subset S$, the residual $\mathbf{r}^{(m)}$ defined in Eq.
(30) can be written in the following form

$\displaystyle\mathbf{r}^{(m)}$ $\displaystyle=$
$\displaystyle\mathbf{y}^{(m)}-\mathbf{X}^{(m)}\hat{\bm{\theta}}^{(m)}=\mathbf{y}^{(m)}-\mathbf{X}_{\mid\hat{S}}^{(m)}\hat{\bm{\theta}}_{\mid\hat{S}}^{(m)}=\left(\mathbf{I}-\mathbf{X}_{\mid\hat{S}}^{(m)}\left(\mathbf{X}_{\mid\hat{S}}^{(m)}\right)^{\dagger}\right)\mathbf{y}^{(m)}$
(78) $\displaystyle=$
$\displaystyle\left(\mathbf{I}-\mathbf{P}_{\hat{S}}^{(m)}\right)\mathbf{y}^{(m)}=\left(\mathbf{I}-\mathbf{P}_{\hat{S}}^{(m)}\right)\left(\tilde{\mathbf{X}}^{(m)}\tilde{\bm{\theta}}^{(m)}+\sigma\bm{\xi}^{(m)}\right)$
$\displaystyle=$
$\displaystyle\left(\mathbf{I}-\mathbf{P}_{\hat{S}}^{(m)}\right)\sum_{l\in
S\setminus\hat{S}}\tilde{\theta}_{l}^{(m)}\tilde{\mathbf{x}}_{l}^{(m)}+\sigma\tilde{\bm{\xi}}^{(m)},$

where $\tilde{\mathbf{X}}^{(m)}$ and $\tilde{\bm{\theta}}^{(m)}$ are the
scaled versions of $\mathbf{X}^{(m)}$ and $\bm{\theta}$, as discussed after
Eq. (8), and the last equality follows from the definition of
$\mathbf{P}_{\hat{S}}^{(m)}$ as a projection operator, so that
$\left(\mathbf{I}-\mathbf{P}_{\hat{S}}^{(m)}\right)\tilde{\mathbf{x}}_{k}^{(m)}=\mathbf{0}$
for any $k\in\hat{S}$.

Denote by $K_{d}$ the size of the detected support set, i.e.,
$K_{d}=\left|\hat{S}\right|$, and by $K_{u}$ the size of the undetected
support set, i.e., $K_{u}=|S\setminus\hat{S}|$. Since $\hat{S}\subset S$, then
$K_{d}+K_{u}=K$. Finally, we introduce the following quantity

$\mu_{\hat{S}}=\frac{K_{d}\mu_{\max}^{2}}{1-(K_{d}-1)\mu_{\max}}.$ (79)

The following Lemma 13 bounds the effect of the projection
$\mathbf{I}-\mathbf{P}_{\hat{S}}^{(m)}$ on the inner products and norms of
columns of $\tilde{\mathbf{X}}^{(m)}$. Its proof appear in Appendix A.6.

######  Lemma 13.

Assume that the max-MIP condition (6) holds and that $\hat{S}\subset S$. Then,
the following inequalities hold:

  1. 1.

The quantity $\mu_{\hat{S}}$ of Eq. (79) satisfies

$\mu_{\hat{S}}\leq\mu_{\max},$ (80)

and

$K_{u}\left(\mu_{\max}+\mu_{\hat{S}}\right)<K\mu_{\max}.$ (81)

  2. 2.

For each index $i\notin\hat{S}$,

$1-\mu_{\hat{S}}\leq\left\|\left(\mathbf{I}-\mathbf{P}_{\hat{S}}^{(m)}\right)\tilde{\mathbf{x}}_{i}^{(m)}\right\|_{2}^{2}\leq
1.$ (82)

  3. 3.

For each pair of distinct indices $i\neq k$ such that $i,k\notin\hat{S}$,

$\left|\left\langle\tilde{\mathbf{x}}_{k}^{(m)},\left(\mathbf{I}-\mathbf{P}_{\hat{S}}^{(m)}\right)\tilde{\mathbf{x}}_{i}^{(m)}\right\rangle\right|\leq\mu_{\max}+\mu_{\hat{S}},$
(83)

and

$\left\|\left(\mathbf{I}-\mathbf{P}_{\hat{S}}^{(m)}\right)\left(\mathbf{I}-\mathbf{P}_{k}^{(m)}\right)\tilde{\mathbf{x}}_{i}^{(m)}\right\|_{2}^{2}\geq
1-\mu_{\max}^{2}-\mu_{\hat{S}}\left(1+\mu_{\max}\right)^{2}.$ (84)

Furthermore, $1-\mu_{\max}^{2}-\mu_{\hat{S}}\left(1+\mu_{\max}\right)^{2}>0.$

For future use, notice that by its definition in Eq. (79), $\mu_{\hat{S}}$ is
an increasing function of $K_{d}$. Since $K_{d}\leq K-1$, then

$\mu_{\hat{S}}\leq\frac{K-1}{1-(K-2)\mu_{\max}}\mu_{\max}^{2}.$ (85)

Thus, under max-MIP condition (6) the quantities $N_{a}$ and $N_{b}$ of Eqs.
(20) and (21) are both strictly positive, and the quantities in Section 4 are
well defined.

For simplicity of notation, from now on we omit the dependence on the machine
index $m$. Given the current estimated support set $\hat{S}$, recall the
definition of $\tilde{\theta}_{\max}$ in Eq. (32) and let $k\in
S\setminus\hat{S}$ be an index for which

$\left\lVert\tilde{\mathbf{x}}_{k}\right\rVert\cdot\left|\theta_{k}\right|=\tilde{\theta}_{\max}$
(86)

(chosen arbitrarily in case of ties). By Eq. (33) for $\rho$ and Eq. (7) for
$\theta_{{\mbox{\tiny crit}}}$,

$\tilde{\theta}_{\max}=\frac{\sigma\sqrt{2\rho\log d}}{1-(2K-1)\mu_{\max}}.$
(87)

We now prove Lemmas 7 and 8.

###### Proof of Lemma 7.

Recall that $q$, defined in Eq. (44), is the probability that some support
index is selected by OMP_Step. A sufficient condition for this to occur is
that the index $k$ defined in Eq. (86) has a higher correlation with the
current residual than any non-support index $j\notin S$. Thus, $q$ is lower
bounded by the probability of the following event

$\left|\left\langle\tilde{\mathbf{x}}_{k},\mathbf{r}\right\rangle\right|\geq\max_{i\notin
S}\left|\left\langle\tilde{\mathbf{x}}_{i},\mathbf{r}\right\rangle\right|.$
(88)

Thus, to prove the lemma it suffices to lower bound the probability of event
(88). Similarly to the proof of Lemma 4, we decompose the noise
$\tilde{\bm{\xi}}$ in Eq. (78) as the sum of two components, the first
$\tilde{\bm{\xi}}_{\parallel}=\mathbf{P}_{k}\tilde{\bm{\xi}}=\left\langle\tilde{\mathbf{x}}_{k},\tilde{\bm{\xi}}\right\rangle\tilde{\mathbf{x}}_{k}$
is parallel to $\tilde{\mathbf{x}}_{k}$, namely
$\left\langle\tilde{\mathbf{x}}_{k},\tilde{\bm{\xi}}_{\parallel}\right\rangle=\left\langle\tilde{\mathbf{x}}_{k},\tilde{\bm{\xi}}\right\rangle$,
and the second
$\tilde{\bm{\xi}}_{\perp}=\tilde{\bm{\xi}}-\tilde{\bm{\xi}}_{\parallel}=\left(\mathbf{I}-\mathbf{P}_{k}\right)\tilde{\bm{\xi}}$,
is orthogonal to $\tilde{\mathbf{x}}_{k}$, i.e.,
$\left\langle\tilde{\mathbf{x}}_{k},\tilde{\bm{\xi}}_{\perp}\right\rangle=0$.

Next, we use this decomposition to bound each of the terms in (88). By Eq.
(78), for any index $i$,

$\left\langle\tilde{\mathbf{x}}_{i},\mathbf{r}\right\rangle=\sum_{l\in
S\setminus\hat{S}}\tilde{\theta}_{l}\left\langle\tilde{\mathbf{x}}_{i},\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\tilde{\mathbf{x}}_{l}\right\rangle+\sigma\left\langle\tilde{\mathbf{x}}_{i},\tilde{\bm{\xi}}\right\rangle.$

For the index $k$,
$\left\|\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\tilde{\mathbf{x}}_{k}\right\|_{2}^{2}=\left\langle\tilde{\mathbf{x}}_{k},\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\tilde{\mathbf{x}}_{k}\right\rangle\geq
1-\mu_{\hat{S}}$ by Eq. (82). For any other undetected support index $l\in
S\setminus\left\\{\hat{S}\cup k\right\\}$,
$\left|\left\langle\tilde{\mathbf{x}}_{k},\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\tilde{\mathbf{x}}_{l}\right\rangle\right|\leq\mu_{\max}+\mu_{\hat{S}}$
by (83) and $\left|\tilde{\theta}_{l}\right|\leq\tilde{\theta}_{\max}$ by its
definition in Eq. (32). Combining these bounds with the definition of
$\tilde{\bm{\xi}}_{\parallel}$ implies that the LHS of (88) can be bounded by

$\displaystyle\left|\left\langle\tilde{\mathbf{x}}_{k},\mathbf{r}\right\rangle\right|$
$\displaystyle\geq$
$\displaystyle\operatorname{sign}\left(\tilde{\theta}_{k}\right)\left\langle\tilde{\mathbf{x}}_{k},\mathbf{r}\right\rangle$
(89) $\displaystyle\geq$
$\displaystyle\tilde{\theta}_{\max}\left(\left\langle\tilde{\mathbf{x}}_{k},\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\tilde{\mathbf{x}}_{k}\right\rangle-\sum_{l\in
S\setminus\left(\hat{S}\cup\left\\{k\right\\}\right)}\left|\left\langle\tilde{\mathbf{x}}_{k},\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\tilde{\mathbf{x}}_{l}\right\rangle\right|\right)+\operatorname{sign}\left(\tilde{\theta}_{k}\right)\sigma\left\langle\tilde{\mathbf{x}}_{k},\tilde{\bm{\xi}}\right\rangle$
$\displaystyle\geq$
$\displaystyle\tilde{\theta}_{\max}\left(1-\mu_{\hat{S}}-\left(K_{u}-1\right)\left(\mu_{\max}+\mu_{\hat{S}}\right)\right)+\sigma\operatorname{sign}\left(\tilde{\theta}_{k}\right)\left\langle\tilde{\mathbf{x}}_{k},\tilde{\bm{\xi}}_{\parallel}\right\rangle.$

The RHS of (88) can be bounded by

$\displaystyle\max_{i\notin
S}\left|\left\langle\tilde{\mathbf{x}}_{i},\mathbf{r}\right\rangle\right|$
$\displaystyle=$ $\displaystyle\max_{i\notin S}\left|\sum_{l\in
S\setminus\hat{S}}\tilde{\theta}_{l}\left\langle\tilde{\mathbf{x}}_{i},\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\tilde{\mathbf{x}}_{l}\right\rangle+\sigma\left\langle\tilde{\mathbf{x}}_{i},\tilde{\bm{\xi}}_{\perp}+\tilde{\bm{\xi}}_{\parallel}\right\rangle\right|$
(90) $\displaystyle\leq$ $\displaystyle\tilde{\theta}_{\max}\max_{i\notin
S}\sum_{l\in
S\setminus\hat{S}}\left|\left\langle\tilde{\mathbf{x}}_{i},\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\tilde{\mathbf{x}}_{l}\right\rangle\right|+\sigma\max_{i\notin
S}\left|\left\langle\tilde{\mathbf{x}}_{i},\tilde{\bm{\xi}}_{\perp}\right\rangle\right|+\sigma\left|\left\langle\tilde{\mathbf{x}}_{k},\tilde{\bm{\xi}}\right\rangle\right|\max_{i\notin
S}\left|\left\langle\tilde{\mathbf{x}}_{i},\tilde{\mathbf{x}}_{k}\right\rangle\right|$
$\displaystyle\leq$ $\displaystyle
K_{u}\tilde{\theta}_{\max}\left(\mu_{\max}+\mu_{\hat{S}}\right)+\sigma\max_{i\notin
S}\left|\left\langle\tilde{\mathbf{x}}_{i},\tilde{\bm{\xi}}_{\perp}\right\rangle\right|+\sigma\mu_{\max}\left|\left\langle\tilde{\mathbf{x}}_{k},\tilde{\bm{\xi}}_{\parallel}\right\rangle\right|.$

where the first step follows from Eq. (78) and the definitions of
$\tilde{\bm{\xi}}_{\perp}$ and $\tilde{\bm{\xi}}_{\parallel}$, the second step
follows from the triangle inequality and the definitions of
$\tilde{\bm{\xi}}_{\parallel}$ and $\tilde{\theta}_{\max}$, and the last
inequality follows from Eq. (83) and the definitions of $\mu_{\max}$ in Eq.
(5). Combining Eq. (89) with Eq. (90) implies that a sufficient condition for
(88) to occur is

$\max_{i\notin
S}\left|\left\langle\tilde{\mathbf{x}}_{i},\tilde{\bm{\xi}}_{\perp}\right\rangle\right|\leq\operatorname{sign}\left(\tilde{\theta}_{k}\right)\left\langle\tilde{\mathbf{x}}_{k},\tilde{\bm{\xi}}_{\parallel}\right\rangle-\mu_{\max}\left|\left\langle\tilde{\mathbf{x}}_{k},\tilde{\bm{\xi}}_{\parallel}\right\rangle\right|+\frac{\tilde{\theta}_{\max}\left(1-2K_{u}\left(\mu_{\max}+\mu_{\hat{S}}\right)+\mu_{\max}\right)}{\sigma}.$

By Eq. (87) and by the inequality (81), a sufficient condition for the
previous event to occur is

$\max_{i\notin
S}\left|\left\langle\tilde{\mathbf{x}}_{i},\tilde{\bm{\xi}}_{\perp}\right\rangle\right|\leq\operatorname{sign}\left(\tilde{\theta}_{k}\right)\left\langle\tilde{\mathbf{x}}_{k},\tilde{\bm{\xi}}_{\parallel}\right\rangle-\mu_{\max}\left|\left\langle\tilde{\mathbf{x}}_{k},\tilde{\bm{\xi}}_{\parallel}\right\rangle\right|+\sqrt{2\rho\log
d}.$ (91)

As in the proof of Lemma 4, denote the LHS of (91) by $A$, its RHS by $B$ and
let $T=\sqrt{2\log d}$. By Eq. (59), it suffices to bound the probabilities of
$A\leq T$ and $B\geq T$.

We begin with bounding $\Pr\left[A\leq T\right]$. Fix $i\notin S$. By
definition
$\tilde{\bm{\xi}}_{\perp}=\left(\mathbf{I}-\mathbf{P}_{k}\right)\tilde{\bm{\xi}}=\left(\mathbf{I}-\mathbf{P}_{k}\right)\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\bm{\xi}$.
By the symmetry of projections,
$\left\langle\tilde{\mathbf{x}}_{i},\tilde{\bm{\xi}}_{\perp}\right\rangle=\left\langle\tilde{\mathbf{x}}_{i},\left(\mathbf{I}-\mathbf{P}_{k}\right)\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\bm{\xi}\right\rangle=\left\langle\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\left(\mathbf{I}-\mathbf{P}_{k}\right)\tilde{\mathbf{x}}_{i},\bm{\xi}\right\rangle.$
Normalizing the inner product results in a standard normal random variable
$Z_{i}=\frac{\left\langle\tilde{\mathbf{x}}_{i},\tilde{\bm{\xi}}_{\perp}\right\rangle}{\left\|\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\left(\mathbf{I}-\mathbf{P}_{k}\right)\tilde{\mathbf{x}}_{i}\right\|_{2}}\sim\mathcal{N}(0,1).$
By Eq. (84),
$\left\|\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\left(\mathbf{I}-\mathbf{P}_{k}\right)\tilde{\mathbf{x}}_{i}\right\|_{2}\geq\gamma_{1},$
where
$\gamma_{1}=\sqrt{1-\mu_{\max}^{2}-\mu_{\hat{S}}\left(1+\mu_{\max}\right)^{2}}$.
As in the proof of Lemma 4, it follows that

$\Pr\left[A\leq T\right]\geq 1-\frac{\gamma_{1}}{\sqrt{\pi\log
d}}d^{1-\frac{1}{\gamma_{1}^{2}}}.$ (92)

We now bound $\Pr\left[B\geq T\right]$, where $B$ is the RHS of Eq. (91). By
definition of $\tilde{\bm{\xi}}_{\parallel}$, the inner product
$\left\langle\tilde{\mathbf{x}}_{k},\tilde{\bm{\xi}}_{\parallel}\right\rangle=\left\langle\tilde{\mathbf{x}}_{k},\tilde{\bm{\xi}}\right\rangle=\left\langle\tilde{\mathbf{x}}_{k},\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\bm{\xi}\right\rangle$.
This random variable is equal in distribution to a Gaussian random variable
$Z\sim\mathcal{N}\left(0,\left\|\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\tilde{\mathbf{x}}_{k}\right\|_{2}^{2}\right)$.
By Eq. (82),
$\left\|\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\tilde{\mathbf{x}}_{k}\right\|_{2}\geq\sqrt{1-\mu_{\hat{S}}}$.
As in the proof of Lemma 4,

$\Pr\left[B\geq
T\right]\geq\Phi^{c}\left(\frac{1-\sqrt{\rho}}{\gamma_{2}}\sqrt{2\log
d}\right),$ (93)

where $\gamma_{2}=\sqrt{1-\mu_{\hat{S}}}\left(1-\mu_{\max}\right)$.

Recall the definitions of $N_{a}$ and $N_{b}$ in Eqs. (20) and (21),
respectively. By Eq. (85),
$\gamma_{1}=\sqrt{1-\mu_{\max}^{2}-\mu_{\hat{S}}\left(1+\mu_{\max}\right)^{2}}\geq\sqrt{N_{b}}$
and
$\gamma_{2}=\sqrt{1-\mu_{\hat{S}}}\left(1-\mu_{\max}\right)\geq\sqrt{N_{a}}(1-\mu_{\max})$.
Combining this with the bounds (92) and (93) completes the proof of Lemma 7. ∎

###### Proof of Lemma 8.

Fix a non-support index $j\notin S$. Recall that $p_{j}$, defined in Eq. (39),
is the probability that index $j$ is selected by OMP_Step. This occurs if $j$
has the highest correlation with the current residual, i.e.,

$p_{j}=\Pr\left[\left|\langle\tilde{\mathbf{x}}_{j},\mathbf{r}\rangle\right|\geq\max_{i\in\left[d\right]\setminus\hat{S}}\left|\langle\tilde{\mathbf{x}}_{i},\mathbf{r}\rangle\right|\right].$

Clearly, by taking the maximum over a subset of the indices
$\mathcal{A}\subseteq\left[d\right]\setminus\hat{S}$ that includes $j$, the
probability can only be higher. Namely,

$p_{j}\geq\Pr\left[\left|\langle\tilde{\mathbf{x}}_{j},\mathbf{r}\rangle\right|\geq\max_{i\in{\cal
A}}\left|\langle\tilde{\mathbf{x}}_{i},\mathbf{r}\rangle\right|\right].$ (94)

Here we take $\mathcal{A}$ as the set of all non-support variables plus the
index $k$, i.e., $\mathcal{A}=\left(\left[d\right]\setminus
S\right)\cup\left\\{k\right\\}$, where $k$ is defined in Eq. (86). Next, we
separately upper bound

$\Pr\left[\left|\left\langle\tilde{\mathbf{x}}_{j},\mathbf{r}\right\rangle\right|>\max_{i\notin
S\cup\left\\{j\right\\}}\left|\left\langle\tilde{\mathbf{x}}_{i},\mathbf{r}\right\rangle\right|\right]$
(95)

and

$\Pr\left[\left|\left\langle\tilde{\mathbf{x}}_{j},\mathbf{r}\right\rangle\right|>\left|\left\langle\tilde{\mathbf{x}}_{k},\mathbf{r}\right\rangle\right|\right]$
(96)

and then upper bound $p_{j}$ using (65) with
$A=\left|\left\langle\tilde{\mathbf{x}}_{j},\mathbf{r}\right\rangle\right|$,
$B=\max_{i\notin\left(S\cup\left\\{j\right)\right\\}}\left|\left\langle\tilde{\mathbf{x}}_{i},\mathbf{r}\right\rangle\right|$,
and
$C=\left|\left\langle\tilde{\mathbf{x}}_{k},\mathbf{r}\right\rangle\right|$.

For later use in both bounds, the random variable $A$ can be upper bounded as
follows

$\displaystyle\left|\left\langle\tilde{\mathbf{x}}_{j},\mathbf{r}\right\rangle\right|$
$\displaystyle=$ $\displaystyle\left|\sum_{l\in
S\setminus\hat{S}}\tilde{\theta}_{l}\left\langle\tilde{\mathbf{x}}_{j},\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\tilde{\mathbf{x}}_{l}\right\rangle+\sigma\left\langle\tilde{\mathbf{x}}_{j},\tilde{\bm{\xi}}\right\rangle\right|$
(97) $\displaystyle\leq$ $\displaystyle\tilde{\theta}_{\max}\sum_{l\in
S\setminus\hat{S}}\left|\left\langle\tilde{\mathbf{x}}_{j},\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\tilde{\mathbf{x}}_{l}\right\rangle\right|+\sigma\left|\left\langle\tilde{\mathbf{x}}_{j},\tilde{\bm{\xi}}\right\rangle\right|\leq\tilde{\theta}_{\max}K_{u}\left(\mu_{\max}+\mu_{\hat{S}}\right)+\sigma\left|\left\langle\tilde{\mathbf{x}}_{j},\tilde{\bm{\xi}}\right\rangle\right|,$

where the first equality follows from Eq. (78), the next inequality follows
from the triangle inequality and the definition of $\tilde{\theta}_{\max}$ in
Eq. (32), and the last inequality follow from (83). We now begin with event
(95). By Eqs. (83) and (81), for each non-support index $i\notin S$ such that
$i\neq j$,

$\displaystyle\left|\left\langle\tilde{\mathbf{x}}_{i},\mathbf{r}\right\rangle\right|$
$\displaystyle\geq$
$\displaystyle\left\langle\tilde{\mathbf{x}}_{i},\mathbf{r}\right\rangle=\sum_{l\in
S\setminus\hat{S}}\tilde{\theta}_{l}\left\langle\tilde{\mathbf{x}}_{i},\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\tilde{\mathbf{x}}_{l}\right\rangle+\sigma\left\langle\tilde{\mathbf{x}}_{i},\tilde{\bm{\xi}}\right\rangle$
$\displaystyle\geq$
$\displaystyle-\tilde{\theta}_{\max}K_{u}\left(\mu_{\max}+\mu_{\hat{S}}\right)+\sigma\left\langle\tilde{\mathbf{x}}_{i},\tilde{\bm{\xi}}\right\rangle.$

Combining the above bound with Eq. (97), rearranging the terms, recalling the
relation between $\tilde{\theta}_{\max}$ and $\rho$ in (87) and applying
inequality (81) yields

$\Pr\left[\left|\left\langle\tilde{\mathbf{x}}_{j},\mathbf{r}\right\rangle\right|>\max_{i\notin
S\cup\left\\{j\right\\}}\left|\left\langle\tilde{\mathbf{x}}_{i},\mathbf{r}\right\rangle\right|\right]\leq\Pr\left[\left|\left\langle\tilde{\mathbf{x}}_{j},\tilde{\bm{\xi}}\right\rangle\right|+\frac{2K\mu_{\max}\sqrt{2\rho\log
d}}{1-(2K-1)\mu_{\max}}>\max_{i\notin
S\cup\left\\{j\right\\}}\left\langle\tilde{\mathbf{x}}_{i},\tilde{\bm{\xi}}\right\rangle\right].$
(98)

As in the proof of Lemma 5, applying (68) with
$T=(1-\epsilon)\sqrt{2\left(1-\mu_{\max}\right)\log d}$ and
$\epsilon\in\left(0,1\right)$ as in (22), we can upper bound (98) by

$\Pr\left[\left|\langle\tilde{\mathbf{x}}_{j},\tilde{\bm{\xi}}\rangle\right|\geq
a\sqrt{2\log d}\right]+\Pr\left[\max_{i\notin
S\cup\left\\{j\right\\}}\langle\tilde{\mathbf{x}}_{i},\tilde{\bm{\xi}}\rangle<(1-\epsilon)\sqrt{2\left(1-\mu_{\max}\right)\log
d}\right],$ (99)

where

$a=(1-\epsilon)\sqrt{1-\mu_{\max}}-\frac{2K\mu_{\max}\sqrt{\rho}}{1-(2K-1)\mu_{\max}}.$

By the symmetry of projection matrices,
$\langle\tilde{\mathbf{x}}_{j},\tilde{\bm{\xi}}\rangle=\langle\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\tilde{\mathbf{x}}_{j},\bm{\xi}\rangle$.
By Eq. (82), the norm
$\left\|\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\tilde{\mathbf{x}}_{j}\right\|_{2}\leq
1$ and thus the first term in (99) is bounded by

$2\Phi^{c}\left(a\sqrt{2\log d}\right).$ (100)

We now bound the second term in (99) using Lemma 12 with
$Z_{i}=\left\langle\frac{\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\tilde{\mathbf{x}}_{i}}{\left\|\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\tilde{\mathbf{x}}_{i}\right\|_{2}},\bm{\xi}\right\rangle$.
Towards this goal, notice that by Eq. (82),
$\left\|\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\tilde{\mathbf{x}}_{i}\right\|_{2}\geq\sqrt{1-\mu_{\hat{S}}}$.
Thus, the second term in (99) is upper bounded by

$\Pr\left[\max_{i\notin
S\cup\left\\{j\right\\}}Z_{i}<\frac{(1-\epsilon)\sqrt{2\left(1-\mu_{\max}\right)\log
d}}{\sqrt{1-\mu_{\hat{S}}}}\right].$ (101)

Furthermore, by Eqs. (82) and (83), for each $i,l\notin
S\cup\left\\{j\right\\}$ such that $i\neq l$,
$\textup{E}\left[Z_{i}Z_{l}\right]\leq\frac{\mu_{\max}+\mu_{\hat{S}}}{1-\mu_{\hat{S}}}.$
Thus, we can apply Lemma 12 with
$\eta=\frac{\mu_{\max}+\mu_{\hat{S}}}{1-\mu_{\hat{S}}}$ and
$\delta=1-\epsilon$ to obtain that (101) is bounded by

$C\left(d-K-1\right)^{-\frac{1-\mu_{\max}}{\mu_{\max}+\mu_{\hat{S}}}\epsilon^{2}}\log^{\left(\frac{1-\mu_{\max}}{\mu_{\max}+\mu_{\hat{S}}}\epsilon-1\right)/2}\left(d-K-1\right).$
(102)

Similarly to the proof of Lemma 5, under condition (22) on $\epsilon$ and for
sufficiently large $d=d(\epsilon)$, (100) is larger than (102), and thus

$\Pr\left[\left|\left\langle\tilde{\mathbf{x}}_{j},\mathbf{r}\right\rangle\right|>\max_{i\notin
S\cup\left\\{j\right\\}}\left|\left\langle\tilde{\mathbf{x}}_{i},\mathbf{r}\right\rangle\right|\right]\leq
4\Phi^{c}\left(a\sqrt{2\log d}\right).$ (103)

We now turn to bounding (96). For the support index $k$, similarly to (89),

$\left|\left\langle\tilde{\mathbf{x}}_{k},\mathbf{r}\right\rangle\right|\geq\operatorname{sign}\left(\tilde{\theta}_{k}\right)\left\langle\tilde{\mathbf{x}}_{k},\mathbf{r}\right\rangle\geq\tilde{\theta}_{\max}\left(1-\mu_{\hat{S}}\right)-\tilde{\theta}_{\max}(K_{u}-1)\left(\mu_{\max}+\mu_{\hat{S}}\right)+\sigma\operatorname{sign}\left(\tilde{\theta}_{k}\right)\left\langle\tilde{\mathbf{x}}_{k},\tilde{\bm{\xi}}\right\rangle.$

Combining the bound above with Eq. (97), recalling the relation between
$\tilde{\theta}_{\max}$ and $\rho$ in (87) and applying the inequality (81)
yields

$\Pr\left[\left|\left\langle\tilde{\mathbf{x}}_{j},\mathbf{r}\right\rangle\right|>\left|\left\langle\tilde{\mathbf{x}}_{k},\mathbf{r}\right\rangle\right|\right]\leq\Pr\left[\left|\left\langle\tilde{\mathbf{x}}_{j},\tilde{\bm{\xi}}\right\rangle\right|-\operatorname{sign}\left(\tilde{\theta}_{k}\right)\left\langle\tilde{\mathbf{x}}_{k},\tilde{\bm{\xi}}\right\rangle>\sqrt{2\rho\log
d}\right].$

We now upper bound this probability. Let
$H=\left\langle\tilde{\mathbf{x}}_{j},\tilde{\bm{\xi}}\right\rangle=\left\langle\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\tilde{\mathbf{x}}_{j},\bm{\xi}\right\rangle$
and
$G=\operatorname{sign}\left(\tilde{\theta}_{k}\right)\left\langle\tilde{\mathbf{x}}_{k},\tilde{\bm{\xi}}\right\rangle=\operatorname{sign}\left(\tilde{\theta}_{k}\right)\left\langle\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\tilde{\mathbf{x}}_{k},\bm{\xi}\right\rangle$.
Notice that by Eqs. (82), (83), and (80), $\sigma_{H}^{2},\sigma_{G}^{2}\leq
1$ and $\left|\sigma_{HG}\right|\leq\mu_{\max}+\mu_{\hat{S}}\leq
2\mu_{\max}$.222Notice that we use the bound $\mu_{\hat{S}}\leq\mu_{\max}$ of
Eq. (80) for simplicity, but in fact combining Eqs. (20) and (85) implies the
tighter bound $\mu_{\hat{S}}\leq 1-N_{a}$, which depends on $K$. Thus, as in
the proof of Lemma 5,

$\Pr\left[\left|H\right|-G>c\right]\leq 2\Phi^{c}\left(\sqrt{\frac{2\rho\log
d}{2+4\mu_{\max}}}\right).$ (104)

By Eq. (65), the probability (94) is at most the minimum between (103) and
(104). By the monotonicity of the Gaussian CDF, (94) is upper bounded by

$4\Phi^{c}\left(\max\left\\{\sqrt{\frac{\rho}{2+4\mu_{\max}}},\left((1-\epsilon)\sqrt{1-\mu_{\max}}-\frac{2K\mu_{\max}\sqrt{\rho}}{1-(2K-1)\mu_{\max}}\right)\right\\}\sqrt{2\log
d}\right).$ (105)

Note that inserting the definitions of $Q_{1}$ and $Q_{2}$ in Eqs. (26) and
(27) respectively, into Eq. (28) and rearranging various terms yields

$\frac{1-\sqrt{r}}{\sqrt{N_{a}}\left(1-\mu_{\max}\right)}+\sqrt{Q_{0}}<\max\left\\{\sqrt{\frac{r}{2+4\mu_{\max}}},\left((1-\epsilon)\sqrt{1-\mu_{\max}}-\frac{2K\mu_{\max}\sqrt{r}}{1-(2K-1)\mu_{\max}}\right)\right\\}.$
(106)

The definitions of $r$ and $\rho$ in Eqs. (10) and (33) imply that $\rho\geq
r$. Thus, $\rho$ satisfies Eq. (28) and hence Eq. (106). The RHS of Eq. (106)
is the same as the maximum in Eq. (105) above. Thus, Eq. (105) is upper
bounded by

$4\Phi^{c}\left(\left(\frac{1-\sqrt{\rho}}{\sqrt{N_{a}}\left(1-\mu_{\max}\right)}+\sqrt{Q_{0}}\right)\sqrt{2\log
d}\right).$ (107)

By the assumption $\rho\leq 1$, we can apply Lemma 11. Hence, by the
definition of $Q_{0}$ in Eq. (25), and by the definition of $F$ in Eq. (23),

$\displaystyle p_{j}$ $\displaystyle\leq$ $\displaystyle
4\Phi^{c}\left(\left(\frac{1-\sqrt{\rho}}{\sqrt{N_{a}}\left(1-\mu_{\max}\right)}+\sqrt{Q_{0}}\right)\sqrt{2\log
d}\right)\leq
4\sqrt{2}d^{-Q_{0}}\Phi^{c}\left(\frac{1-\sqrt{\rho}}{\sqrt{N_{a}}\left(1-\mu_{\max}\right)}\sqrt{2\log
d}\right)$ $\displaystyle=$ $\displaystyle
4\sqrt{2}\frac{1-\sqrt{\frac{N_{b}}{\pi\log
d}}d^{1-1/N_{b}}}{44\sqrt{2}K}\Phi^{c}\left(\frac{1-\sqrt{\rho}}{\sqrt{N_{a}}\left(1-\mu_{\max}\right)}\sqrt{2\log
d}\right)=\frac{F\left(d,K,\mu_{\max},\rho\right)}{11K},$

which completes the proof of Lemma 8. ∎

###  A.6 Proofs of technical lemmas

###### Proof of Lemma 11.

Classical results by Birnbaum [9] and Komatu [25] are that for all $x\geq 0$,
the following inequalities hold

$\frac{2e^{-x^{2}/2}}{\sqrt{2\pi}\left(\sqrt{x^{2}+4}+x\right)}<\Phi^{c}\left(x\right)<\frac{2e^{-x^{2}/2}}{\sqrt{2\pi}\left(\sqrt{x^{2}+2}+x\right)}.$
(108)

Hence,

$\Phi^{c}\left(a+b\right)<\frac{2e^{-(a+b)^{2}/2}}{\sqrt{2\pi}\left(\sqrt{(a+b)^{2}+2}+a+b\right)},$

and

$\Phi^{c}\left(a\right)>\frac{2e^{-a^{2}/2}}{\sqrt{2\pi}\left(\sqrt{a^{2}+4}+a\right)}.$

Combining the two yields the following

$\Phi^{c}\left(a+b\right)<\frac{\left(\sqrt{a^{2}+4}+a\right)e^{-ab}}{\left(\sqrt{(a+b)^{2}+2}+a+b\right)}e^{-b^{2}/2}\Phi^{c}\left(a\right).$

Notice that for any $a\geq 0$, the fraction in the above display is a
decreasing function of $b$. Since $b\geq 0$, it suffices to note that

$\frac{\left(\sqrt{a^{2}+4}+a\right)}{\left(\sqrt{a^{2}+2}+a\right)}\leq\sqrt{2}$
for any $a\geq 0$. ∎

Towards proving Lemma 13, we prove the following Lemma 14, which bounds the
inner product between vectors projected to the subspace orthogonal to
$\tilde{\mathbf{X}}_{\mid\hat{S}}$ under the assumption $\hat{S}\subset S$. To
this end, recall that we denote by $K_{d}$ the number of detected support
indices, i.e., $K_{d}=|\hat{S}|$, by $K_{u}$ the number of undetected support
indices, i.e., $K_{u}=|S\setminus\hat{S}|$. Since $\hat{S}\subset S$, then
$K_{d}+K_{u}=K$ and there is at least one undetected support index, hence
$1\leq K_{u}\leq K$ and $0\leq K_{d}\leq K-1$.

######  Lemma 14.

If $\hat{S}\subset S$ and $\left(K_{d}-1\right)\mu_{\max}<1$, then for any
pair of vectors $\mathbf{a}_{1},\mathbf{a}_{2}\in\mathbb{R}^{n}$

$\left\langle\mathbf{a}_{1},\mathbf{a}_{2}\right\rangle-\frac{\left|\sum_{j\in\hat{S}}\left\langle\tilde{\mathbf{x}}_{j},\mathbf{a}_{1}\right\rangle\left\langle\tilde{\mathbf{x}}_{j},\mathbf{a}_{2}\right\rangle\right|}{1-\left(K_{d}-1\right)\mu_{\max}}\leq\left\langle\mathbf{a}_{1},\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\mathbf{a}_{2}\right\rangle\leq\left\langle\mathbf{a}_{1},\mathbf{a}_{2}\right\rangle+\frac{\left|\sum_{j\in\hat{S}}\left\langle\tilde{\mathbf{x}}_{j},\mathbf{a}_{1}\right\rangle\left\langle\tilde{\mathbf{x}}_{j},\mathbf{a}_{2}\right\rangle\right|}{1-\left(K_{d}-1\right)\mu_{\max}}.$
(109)

If in addition $\mathbf{a}_{1}=\mathbf{a}_{2}=\mathbf{a}$, then

$\left\|\mathbf{a}\right\|_{2}^{2}-\frac{\sum_{j\in\hat{S}}\left\langle\tilde{\mathbf{x}}_{j},\mathbf{a}\right\rangle^{2}}{1-\left(K_{d}-1\right)\mu_{\max}}\leq\left\|\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\mathbf{a}\right\|_{2}^{2}\leq\left\|\mathbf{a}\right\|_{2}^{2}-\frac{\sum_{j\in\hat{S}}\left\langle\tilde{\mathbf{x}}_{j},\mathbf{a}\right\rangle^{2}}{1+\left(K_{d}-1\right)\mu_{\max}}.$
(110)

###### Proof of Lemma 14.

First, if $\hat{S}=\emptyset$ then clearly
$\left\langle\mathbf{a}_{1},\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\mathbf{a}_{2}\right\rangle=\left\langle\mathbf{a}_{1},\mathbf{a}_{2}\right\rangle$
and both (109) and (110) trivially hold. Therefore, assume that $K_{d}\geq 1$.
In this case

$\left\langle\mathbf{a}_{1},\left(\mathbf{I}-\mathbf{P}_{\hat{S}}\right)\mathbf{a}_{2}\right\rangle=\left\langle\mathbf{a}_{1},\mathbf{a}_{2}\right\rangle-\left\langle\mathbf{a}_{1},\mathbf{P}_{\hat{S}}\mathbf{a}_{2}\right\rangle.$
(111)

By definition of $\mathbf{P}_{\hat{S}}$ in Eq. (50),

$\left\langle\mathbf{a}_{1},\mathbf{P}_{\hat{S}}\mathbf{a}_{2}\right\rangle=\left\langle\mathbf{a}_{1},\tilde{\mathbf{X}}_{\mid\hat{S}}\left(\tilde{\mathbf{X}}_{\mid\hat{S}}^{T}\tilde{\mathbf{X}}_{\mid\hat{S}}\right)^{-1}\tilde{\mathbf{X}}_{\mid\hat{S}}^{T}\mathbf{a}_{2}\right\rangle=\left\langle\tilde{\mathbf{X}}_{\mid\hat{S}}^{T}\mathbf{a}_{1},\left(\tilde{\mathbf{X}}_{\mid\hat{S}}^{T}\tilde{\mathbf{X}}_{\mid\hat{S}}\right)^{-1}\tilde{\mathbf{X}}_{\mid\hat{S}}^{T}\mathbf{a}_{2}\right\rangle.$
(112)

We now bound this term in absolute value. For a matrix
$\mathbf{A}\in\mathbb{R}^{n\times n}$, denote by
$\lambda_{\min}\left(\mathbf{A}\right)$ and
$\lambda_{\max}\left(\mathbf{A}\right)$ its minimal and maximal eigenvalues,
respectively. Consider