# SPriFed-OMP: A Differentially Private Federated Learning Algorithm for
Sparse Basis Recovery

Ajinkya K Mulay<EMAIL_ADDRESS>
Elmore Family School of Electrical and Computer Engineering
Purdue University Xiaojun Lin<EMAIL_ADDRESS>
Elmore Family School of Electrical and Computer Engineering
Purdue University

###### Abstract

Sparse basis recovery is a classical and important statistical learning
problem when the number of model dimensions $p$ is much larger than the number
of samples $n$. However, there has been little work that studies sparse basis
recovery in the Federated Learning (FL) setting, where the client data’s
differential privacy (DP) must also be simultaneously protected. In
particular, the performance guarantees of existing DP-FL algorithms (such as
DP-SGD) will degrade significantly when $p\gg n$, and thus, they will fail to
learn the true underlying sparse model accurately. In this work, we develop a
new differentially private sparse basis recovery algorithm for the FL setting,
called SPriFed-OMP. SPriFed-OMP converts OMP (Orthogonal Matching Pursuit) to
the FL setting. Further, it combines SMPC (secure multi-party computation) and
DP to ensure that only a small amount of noise needs to be added in order to
achieve differential privacy. As a result, SPriFed-OMP can efficiently recover
the true sparse basis for a linear model with only $n=\mathcal{O}(\sqrt{p})$
samples. We further present an enhanced version of our approach, SPriFed-OMP-
GRAD based on gradient privatization, that improves the performance of
SPriFed-OMP. Our theoretical analysis and empirical results demonstrate that
both SPriFed-OMP and SPriFed-OMP-GRAD terminate in a small number of steps,
and they significantly outperform the previous state-of-the-art DP-FL
solutions in terms of the accuracy-privacy trade-off.

##  1 Introduction

For many statistical learning applications, such as genomics and economics
(Liang & Kelemen, 2008; Clarke et al., 2008; Belloni et al., 2014; Fan et al.,
2011), it is essential to deal with situations where the number of samples
($n$) is significantly lower than the number of model parameters ($p$).
Without additional constraints, such problems are ill-determined as there will
be many models that can fit the same set of training samples. To address this
issue, a popular approach is to assume some sparsity conditions on the desired
model. Along this line, there have been significant advances in compressed
sensing techniques, such as LASSO (Zhao & Yu, 2006), (Meinshausen, 2007)
Orthogonal Matching Pursuit (OMP) (Tropp et al., 2007), (Needell & Tropp,
2009) Forward-Backward Algorithm (FoBA) (Zhang, 2011) and Least Angles
Regression (LARS) (Efron et al., 2004), those leverage sparsity assumptions on
the data to extract the true sparse basis of the underlying model. Once the
sparse basis is identified, the sparse model can be estimated using standard
learning techniques, as the problem is no longer ill-determined.

However, the above methods usually assume that the training datasets are
already hosted on a server and thus neglect the privacy of the users who
generate the training data. Recently, Federated Learning (FL) (McMahan et al.,
2017) has been proposed as a standard approach to protect user privacy, which
keeps user data at the clients during training. Specifically, the clients
update the models on their data and only upload the new gradients to the
parameter server. The server updates the model parameters based on the
aggregated global gradient and then sends the new model parameters back to the
clients, furthering the training process. Such iterations are repeated until
the global model converges, with no need to upload the user data to the
server. However, note that FL is insufficient for privacy protection because
the client gradients transmitted to the server may still leak information
(Huang et al., 2021), (Geiping et al., 2020). Therefore, in the literature,
Differential Privacy (DP) (Dwork et al., 2014) has also been introduced into
FL. DP privatizes FL by adding noise to client gradients before uploading them
to the parameter server. By adding appropriate noises, individual client
gradients are hidden from the server and any external adversary while still
enabling the learning of the global model. A notable example is DP-SGD, which
can be applied to any SGD-based training algorithms and arbitrary loss
functions (e.g., regularized objective functions in LASSO).

Unfortunately, this DP-FL pipeline is a poor fit for training sparse models
when the number of training samples ($n$) is significantly below the number of
model parameters ($p$). For ease of discussion, in this paper, we assume that
each client contains a single training sample. Therefore, in our analysis, the
number of clients equals the number of samples. 111Our analysis can be easily
extended to multiple samples per client by the group privacy notion discussed
in Section 2. Note that, intuitively, the gradient uploaded by each client has
$p$ elements, one for each model parameter. For DP to achieve the desirable
privacy guarantee, we usually need to add noise to the gradients with variance
proportional to the model dimensions (Dwork et al., 2014). Thus, when $p\gg
n$, the noise required in the DP-FL setting can completely overwhelm the
signal and thus prevent the recovery of the correct sparse basis. We refer to
this problem as the curse-of-dimensionality for DP-FL.

In particular, for DP-SGD (Abadi et al., 2016), even when it is applied to
Lipschitz loss functions, the empirical risk is of the order
$\mathcal{O}(\frac{p}{n})$ (Theorem 2.4 in Bassily et al. (2014)). When $p\gg
n$, even though one can apply DP-SGD to a sparse-solution-seeking objective,
such as LASSO, it will not produce accurate answers. Other DP-FL algorithms,
such as objective perturbation (Kifer et al., 2012), have similar issues. For
example, the empirical loss of the objective perturbation mechanism in Kifer
et al. (2012) is of the order $\mathcal{O}(\frac{p^{2}}{n})$ (see Theorems
4.1-4.2 and 5 in Kifer et al. (2012)), which will also fail to produce
accurate sparse models when $p\gg n$. Thus, developing DP-FL algorithms that
can attain provable sparse recovery under high-dimensional settings remains an
open question. We note that if we do not consider the FL setting, there are DP
solutions in the literature for sparse recovery under high dimension (Thakurta
& Smith, 2013; Talwar et al., 2015). However, these algorithms assume full
data access at the server and thus fail to work in the FL setting (where data
are kept at the client for privacy).

In this work, we develop a new sparse-recovery algorithm specifically for the
DP-FL setting, which, under suitable conditions, can ensure DP for data at the
clients and recover the true sparse basis even when $p\gg n$. Specifically,
our new algorithm called SPriFed-OMP (Sparse Private Federated OMP), is based
on Orthogonal Matching Pursuit (OMP) (Tropp et al., 2007). At each iteration,
OMP picks one basis with the highest correlation with the target, subtracts
the contribution of this basis from the target, and continues to search for
the next basis with the highest correlation with the residual until a given
number $s$ of sparse basis is identified. However, the standard OMP (Tropp et
al., 2007) is not designed for the FL setting, nor does it ensure DP. We
augment OMP with a noisy SMPC (secure multi-party computation) algorithm
(Bonawitz et al., 2017; Kairouz et al., 2021) so that only a smaller amount of
noise is added to the aggregate correlation, which ensures DP. In some steps,
e.g., to compare the correlation with the target across all bases, this noise
is still in order $O(\sqrt{p})$. However, we only need to perform such steps
$s$ times, which is usually much smaller than the number of iterations for DP-
SGD to converge. As a result, the overall amount of noise will be much
smaller. Our careful analysis shows that, as long as
$n=\mathcal{O}(\sqrt{p})$, our algorithm will be able to recover the true
sparse basis with high probability under the similar assumption of Restricted
Isometry Property (RIP) as in standard OMP. While this general idea is quite
intuitive, the detailed design of SPriFed-OMP also matters a lot. Indeed, we
present two versions of SPriFed-OMP. In the first version, we privatize the
individual correlations (feature and feature-residual correlations) to compute
the total feature-residual correlation in each step. This version introduces
our first enhancement that specifically adds lower order noise to the selected
basis elements. Lower order noise enables lower sample size requirement for
sparse basis recovery while also incurring lower test error as we will see in
Section 7. Next, we identify a simplified second version of SPriFed-OMP that
introduces our second enhancement based on gradient privatization. We will see
in the experimental section (Section 7) that adding noise to the gradient is
much more advantageous under clipping than adding noise to correlations.
Furthermore, for both methods, we can quantify the estimation error (Theorem
9) and the empirical risk (Theorem 10) to be on the order of
$\mathcal{O}\Big{(}\sqrt{\frac{s\log(s)}{n}}\Big{)}$ and
$\mathcal{O}\Big{(}\frac{s\log(s)}{n}\Big{)}$, respectively which are in the
same order as the traditional non-private Ordinary Least Squares (OLS)
estimate that already assumes knowledge of the correct sparse basis.

###  1.1 Related Work

In the literature, we have two major categories for private sparse learning:
(1) The sparse learning methods that only work in the central DP setting
(_i.e.,_ where the data can be sent to the server non-privately and only the
server privatizes the output) and, (2) the learning methods that are either
directly DP-FL compatible or can be adapted to be suitable for the DP-FL
setting. We note that not all of these methods can provide basis recovery
guarantees (some of them only guarantees empirical risks). Below, we discuss
the related work in each category.

Central DP compatible sparse learning methods: The recent sparse methods in
Asi et al. (2021); Bassily et al. (2021) have empirical risk that holds for
$n=\mathcal{O}(\sqrt{p})$ samples under the RSC constraint. However, both of
these methods require centralized data access due to their unique sampling
requirements and certain non-private mechanisms, and thus, they are not
suitable for the DP-FL domain.

DP-FL compatible sparse learning methods: As we discussed earlier, DP-SGD is
not designed for the $p\gg n$ setting or exploiting sparsity. Amongst the DP-
FL sparse learning methods, the work of Kifer et al. (2012) is most closely
related to ours. Although the method proposed in Kifer et al. (2012) is not
geared toward the DP-FL setting, the Samp-Agg algorithm can be adapted to the
DP-FL setting. However, support recovery in Kifer et al. (2012) requires that
each client’s data matrix individually satisfies the Restricted Strong
Convexity (RSC) (Zhang, 2011) assumption. In contrast, the RIP assumption in
our proposed solution is considerably weaker, as it only needs to hold over
the entire training set across all clients. Thus, our solution can apply to a
larger set of scenarios. There are two other related work that aim at the DP-
FL and $p\gg n$ setting. Mangold et al. (2023) recently proposed a mechanism
that has an empirical risk that holds for $n=\mathcal{O}(\sqrt{p})$ samples
under the RSC constraint. We note that Mangold et al. (2023) is only studied
for the centralized DP setting. Still, the algorithm can be appropriately
adjusted for the DP-FL setting with the Gaussian noise 222We provide a
modified/enhanced DP-FL version of Algorithm 1 in Mangold et al. (2023) for
comparison in the form of Algorithm 5.. However, their analysis lacks a sparse
basis recovery guarantee, and thus it is unclear whether their methods can
find the underlying true basis. Furthermore, their empirical risk values
cannot avoid the polynomial dependency on $p$. Wang & Xu (2019) also considers
the problem of estimating sparse models under the $p\gg n$ scenario and where
noise is directly added to the client communication to the server to preserve
DP. They provide two algorithmic results. The first result applies to the
$p\gg n$ case, but it studies the privacy scenarios where only the response
vector $\bm{y}$ must be private. The measurement matrix $\bm{X}$ has no
privacy protection. When the matrix $\bm{X}$ must also be protected, their
second result still requires $n=O(p)$ (ignoring the logarithmic terms). Thus,
our proposed algorithms provide better privacy guarantees with fewer samples.
In a similar sense, the private hard-thresholding methods proposed in Wang &
Gu (2019); Hu et al. (2022) either do not recover the sparse models or require
$n=O(p)$ samples for convergence.

Noisy OMP: Finally, our analysis of the noisy version of OMP is also related
to (Chen & Caramanis, 2012), which studies how adding noise to the measurement
matrix $\bm{X}$ will impact the accuracy of OMP. However, their analysis does
not apply to the DP or the DP-FL cases as it does not privatize the response
vector $\bm{y}$. Furthermore, they directly add noise to the measurement
matrix $\bm{X}$. As a result, the noise level needed for ensuring DP for the
matrix $\bm{X}$ would have had to be at least $\sqrt{pn}$, which is higher
than our solution’s. Specifically, our solution adds noise to the correlation
between each basis and the target. Further, by employing Noisy-SMPC (Algorithm
1), the level of noise added is only of the order $\sqrt{ps}$. Thus, a new
analysis is needed to study our improved DP and accuracy guarantees.

The rest of the paper is structured as follows. Section 2 presents the system
model and provides a brief overview of differential privacy and the goal of
compressed sensing. In Section 3, we provide the details of our proposed
algorithm SPriFed-OMP and its other variant SPriFed-OMP-GRAD. Sections 4 and 5
provide a thorough utility-privacy analysis of our proposed implementation.
Section 6 provides an intuitive sketch of the core result. Section 7 provides
empirical results. Then, we conclude. Complete proofs and supporting results
are omitted due to the page limits.

##  2 System Model

We assume that the ground truth model is
$\bm{y}=\bm{x}\bm{\alpha_{*}}+\bm{\epsilon}$, where
$\bm{\alpha_{*}}\in\mathbb{R}^{p}$ is the underlying model parameter that we
wish to recover, $\bm{x}\in\mathbb{R}^{p}$ and $\bm{y}\in\mathbb{R}$ are the
input-output pair, and
$\bm{\epsilon}\sim\mathcal{N}(0,\sigma^{2}_{\epsilon})\in\mathbb{R}$ is an
additive error. We assume that this ground-truth model is sparse so that
$\bm{\alpha_{*}}$ has at most $s$ non-zero elements, i.e.,
$||\bm{\alpha_{*}}||_{0}\leq s$. We now consider $n$ input-output training
pairs represented by
$(\bm{X},\bm{y})\triangleq(\bm{x_{i}},y_{i})_{i\in\\{1,...,n\\}}$ where each
pair belongs to a single distinct client. 333This assumption can be easily
relaxed to allow more than one training pair per client by performing the
privacy analysis with group privacy. As long as the samples of each client are
considerably lower than $n$, the group privacy cost will be small. If a client
holds a large number of samples, we could consider limiting the maximum number
of samples per client in each iteration of our algorithm to upper bound our
privacy costs. We emphasize that the server does not have the entire data, so
we need to study a federated learning setting (McMahan et al., 2017).

Any DP-FL algorithm iterates over $n$ distributed clients over steps
$t=1,...,T$. For each step $t$, the $i^{th}$ client transmits message
$m_{i,t}$, computed as a function of the client dataset and any information it
receives from the server plus possibly additional randomness. At the end of
the $T$-th iteration, the server aims to recover $\bm{\alpha_{*}}$ (likely
with some error). Let us use $\mathcal{M}$ to denote such a distributed and
randomized mechanism. We aim to ensure Differential Privacy (DP) for the
clients, and hence we present the Differential Privacy (DP) definition for
$\mathcal{M}$ below.

######  Definition 2.1.

[Approximate Distributed Differential Privacy] For any two neighboring
distributed datasets $\bm{X,X{{}^{\prime}}}\in\mathbb{R}^{n\times p}$
differing in one of their clients and the corresponding data-pairs (_i.e._ ,
they only differ in the data-pair for one client $k\in\\{1,...,n\\}$). Let
$\mathcal{M}(\bm{X})$ be the output when $\mathcal{M}$ operates on
$([\bm{x_{1},...,x_{k},...,x_{n}}],[y_{1},...,y_{k},...y_{n}])$ and let
$\mathcal{M}(\bm{X^{{}^{\prime}}})$ be the output when $\mathcal{M}$ operates
on
$([\bm{x_{1},...,x_{k}^{{}^{\prime}},...,x_{n}}],[y_{1},...,y_{k}^{{}^{\prime}},...,y_{n}])$.
We say that the randomized and distributed mechanism $\mathcal{M}$ is
$(\epsilon,\delta)$ differentially private if, for any set
$\mathcal{F}\in\mathbb{R}^{p}$,

$\displaystyle\text{Pr}[\mathcal{M}(\bm{X})\in\mathcal{F}]\leq
e^{\epsilon}\cdot\text{Pr}[\mathcal{M}(\bm{X^{{}^{\prime}}})\in\mathcal{F}]+\delta.$

Next, we introduce essential assumptions relevant to our system model. Recall
that the goal of the server is to estimate the sparse model parameter
$\bm{\alpha_{*}}$ without access to the client data. Even without the DP
requirement, recovering the true sparse model typically requires additional
assumptions on the data matrix (Zhang, 2009; Candes & Tao, 2005). Below, in
assumption 1 we describe the Restricted Isometry Property (RIP), which is
usually required for OMP (Candes & Tao, 2005), which our proposed algorithm is
based on.

######  Assumption 1.

Restricted Isometry Property (RIP): A measurement matrix
$\bm{X}\in\mathbb{R}^{n\times p}$ is said to satisfy RIP of order-$K$ if there
exists a constant $\zeta\in[0,1]$ such that,

$\displaystyle(1-\zeta)||\bm{v}||_{2}^{2}\leq||\bm{Xv}||_{2}^{2}\leq(1+\zeta)||\bm{v}||_{2}^{2}$

for all vectors $\bm{v}\in\mathbb{R}^{p},||\bm{v}||_{0}\leq K$. Note that
$\zeta$ is referred to as the isometric constant. Furthermore, the Restricted
Isometry Constant (RIC) is defined as the infimum of all possible $\zeta$
values that satisfy RIP of order-$K$ for a given measurement matrix
$\bm{X}\in\mathbb{R}^{n\times p}$. In other words, assuming that the RIP of
order $K$ is satisfied, the corresponding RIC is given by,

$\displaystyle\zeta_{K}=\text{inf}_{\zeta}$
$\displaystyle\Big{\\{}\zeta|(1-\zeta)||\bm{v}||_{2}^{2}\leq||\bm{Xv}||_{2}^{2}\leq(1+\zeta)||\bm{v}||_{2}^{2}$
$\displaystyle\text{ for all }\bm{v}\in\mathbb{R}^{p},||\bm{v}||_{0}\leq
K\Big{\\}}.$

We acknowledge that RIP is a more restrictive condition than other related
conditions in the compressed-sensing literature, such as the Restricted Strong
Convexity from Jalali et al. (2011) and Zhang (2009) or the Positive Cone
Condition from Efron et al. (2004). However, even under this RIP assumption,
no solution exists in the literature to attain DP in an FL setting and achieve
exact support recovery with the number of samples much smaller than the model
dimension. Thus, our contribution under the RIP assumption still represents a
significant contribution. We will leave the study of other assumptions for
future work.

We also make the following assumption, typical in the differential privacy
literature (Dwork et al., 2014; Wang & Xu, 2019).

######  Assumption 2.

[Bounded data matrix and response] We assume that the elements in the
measurement matrix $\bm{X}$ and its response $\bm{y}$ are bounded by scalars
$X_{M}$ and $y_{M}$ respectively.

The boundedness of matrix or vector elements can be easily achieved by
clipping the values to particular bounds. We can also easily maintain the
original vector variance by re-scaling after clipping. Alternatively, if the
underlying data distribution is light-tailed, we could leverage concentration
bounds combined with the union bound to obtain a realistic bound on the data
values.

As discussed in section 1, the standard DP approaches of adding noise are
highly ineffective for recovering sparse models when $p\gg n$. Thus, our paper
aims to develop a mechanism that is both DP-FL and can recover the exact
support with several samples much smaller than the total model dimension.
Specifically, let $\bm{\hat{\alpha}}$ denote the estimated model parameter of
our proposed DP-FL algorithm. Recall that $\bm{\alpha_{*}}$ is the ground-
truth model parameter with sparsity $s$. We wish to (1) Ensure that the
support of $\bm{\hat{\alpha}}$ matches the true support with high probability,
(2) Quantify the empirical risk $\Delta
R\triangleq\frac{1}{n}\sum_{i=1}^{n}\Big{(}(\bm{x_{i}}\bm{\hat{\alpha}}-y_{i})^{2}-(\bm{x_{i}}\bm{\alpha_{*}}-y_{i})^{2}\Big{)}$
and estimation error
$\Delta\alpha\triangleq||\bm{\hat{\alpha}}-\bm{\alpha_{*}}||_{2}$ such that
both are small and do not blow up with large $p$.

##  3 The SPriFed-OMP Algorithm

Our primary goal in this work is to recover a sparse model given a high number
of features and a few samples. Note that if we can guarantee exact sparse
basis recovery, we only need to add noise to the output model with variance
proportional to the model sparsity, which is a much easier goal to accomplish.
Thus, below, we will first focus on the goal to identify the true sparse
basis. The analysis for recovering $\alpha_{*}$ with a small estimation error
will then be presented later (see Theorems 9 and 10).

Non-Private OMP: Our proposed algorithm is based on OMP, a popular algorithm
for exact sparse recovery without DP considerations. Below, we first describe
the standard version of OMP. OMP iteratively selects a single new feature in
each step. The feature with the absolute maximum correlation to the current
residual is picked during each step. The residual is set to the model response
$y$ in the initial step. At each subsequent step, the residual is updated by
removing the newly picked feature’s contribution from the response. OMP
continues iterating until a model of a predetermined dimension/sparsity $s$ is
selected.

Challenge to make OMP differentially private: Below, we explain the challenge
to make OMP differentially private in the FL setting. First, when OMP computes
the correlation of each basis with the residual and picks the basis with the
highest correlation, it must be able to do so without direct access to the
data on all clients. Second, when OMP subtracts the contribution of a new
basis from the current residual, it needs the covariance of the existing
basis. This must also be done without direct access to the data on all
clients.

Both of these challenges can be overcome if we can differentially privately
compute the correlation between two columns that are spread across all
clients. Further, the total number of such computations must be carefully
controlled to limit the amount of noise added. This idea leads to our first
version of SPriFed-OMP, which focuses on computing the correlations in a DP-FL
manner. The server maintains the estimated model parameter $\hat{\alpha}$ and
is never released to the clients. (In contrast, the second version of SPriFed-
OMP, which will be presented later, will release the privatized model
parameter to the clients.)

Differentially-Privately and Distributed Computation of Correlation :
Specifically, consider two columns of data $z^{a}_{i}$ and $z^{b}_{i}$ spread
across all clients $i\in\\{1,...,n\\}$. We now develop a mechanism so that the
server can compute the sum $\sum_{i=1}^{n}z^{a}_{i}z^{b}_{i}$. For privacy,
though, each client cannot disclose $z^{a}_{i}z^{b}_{i}$ directly. We can add
noise directly to the client, which will lead to a noise level that is too
high. In contrast, below, we use an approach that requires lower noise.

NoisySMPC: We use the NoisySMPC mechanism (Algorithm 1) that combines SMPC
(Bonawitz et al., 2017) with a much lower level of noise to ensure DP. This
algorithm, which is distributed and satisfies DP, forms a core component of
our proposed SPriFed-OMP algorithm. To achieve DDP (distributed and
differentially private), NoisySMPC adds two levels of randomness. First,
NoisySMPC modifies each client’s contribution as
$\tilde{f}(\bm{x_{i}},y_{i})=f(\bm{x_{i}},y_{i})+\eta_{i}$ where
$\bm{x_{i}},y_{i}$ is the $i^{th}$ client’s data-response pair and $f$
computes a statistic such as the correlation between the client’s data and
response. The per-client noise $\eta_{i}$ is added so that the total noise
$\eta_{sum}=\sum_{i=1}^{n}\eta_{i}$ at the server is sufficient to
differentially privatize the sum of client outputs. Second, NoisySMPC adds
$f(z^{a}_{i},z^{b}_{i})$ across all clients through SMPC, which further
protects the privacy of individual clients. This underlying SMPC mechanism
allows the clients to sum their contributions $\tilde{f}$ without disclosing
any individual values. Examples of such SMPC mechanisms can be pair-wise
client key sharing as described in (Bonawitz et al., 2017) or the distributed
discrete Gaussian (Kairouz et al., 2021). We refer the reader to the above
literature for further details regarding SMPC and related mechanisms.
NoisySMPC adds significantly lower variance noise to the sum (reduced by a
factor of $n$), compared to privatizing each client’s output individually, as
shown in line $11$ of the Algorithm 1).

Private OMP: We are now ready to present the complete SPriFed-OMP algorithm in
the DP-FL setting. Algorithm 3 contains the pseudocode. Recall that, for each
step in OMP, we compute, for all features, the correlation of the feature
column and the residual. For the first step, the residual is set to the model
response, and thus, the computation can be denoted by $\bm{X^{T}y}$. The lines
$4-6$ perform this computation privately using the NoisySMPC Algorithm 1. This
noisy correlation selects the feature with the highest absolute correlation as
part of the predicted basis (line $8$). We then subtract the correlation
contributed by the newly contributed feature $l^{*}$ from the previous
residual to obtain the new residual. Mathematically, we represent this
correlation contribution as
$\bm{\beta_{l^{*}}(\beta_{l^{*},l^{*}})^{-1}\gamma_{l^{*}}}=\bm{X^{T}X_{l^{*}}(X_{l^{*}}^{T}X_{l^{*}})^{-1}X_{l^{*}}^{T}y}$
where
$\bm{\beta_{l^{*}}=X^{T}X_{l^{*}},\beta_{l^{*},l^{*}}=(X_{l^{*}}^{T}X_{l^{*}}),\gamma_{l^{*}}=X_{l^{*}}^{T}y}$.
Note that this computation computes the correlation between the new column
$l^{*}$ and other columns. These correlations can be computed via the
NoisySMPC mechanism. We represent these combined computations on lines
$10-18$. Note that after a feature is chosen, its associated correlation
computation must only be done twice via NoisySMPC. We first privatize
correlations over all un-selected features (order $\sqrt{p-l}$ where $l$ are
the number of features already chosen). From amongst these correlations, we
choose the correlation with the highest value. Once the highest correlated
feature is identified, we add noise to this correlation with a much lower
noise value (of order $\mathcal{s}$ where $s$ is the sparsity of the model).
Thus, we can significantly reduce the noise impact by privatizing the highest
correlation twice.

_Remark:_ We notice that, to privatize the $p$-dimensional correlations (lines
5 and 11), we cannot reduce the variance of DP noise below $\mathcal{O}(p)$
(see the Gaussian mechanism from Dwork et al. (2014)). Thus, each time we
compute the $p$-dimensional correlations, we are required to add noise of
magnitude $\mathcal{O}(\sqrt{p})$. However, we only need to do so finitely
many times (_i.e.,_ $s$ times), thus requiring a finite privacy budget.
Further, thanks to NoisySMPC, we can get away with adding a significantly
smaller amount of noise (lines 14-17).

Further, we note that for our method, we need to compute a maximum of
$p(s+1)+2s^{2}$ private correlations. In DP-SGD, we instead need to compute
$pT$ private correlations, where the number of iterations $T$ is significantly
higher than $s$ and often around the order of $n$. Furthermore, in our
proposed algorithm, we require order-$p$ variance noise only for recovering
the basis while order-$s$ variance noise is used to compute the model
parameter on these bases. DP-SGD on the other hand requires order-$p$ variance
noise throughout for all $T$ iterations. Thus, the noise required for our
method is much lower in terms of the number of private correlation
computations and the model parameter computation. Therefore, we expect to
recover a model with significantly higher accuracy.

_Enhancement 1: Adding lower noise to selected features_

While the above idea of privatize OMP may seem straightforward, there is a key
step in Line 14 which greatly enhances the performance of SPriFed-OMP. Note
that in lines 5 and 11 of Algorithm 3, we already privatize the $p$
dimensional artifacts $\bm{X^{T}y}$ and $\bm{X^{T}X_{l^{*}}}$ respectively
with noise of variance $p\sigma_{1}^{2}$. Although we could have directly use
such privatized values in estimating $\bm{\alpha}$, the resulting error will
be high. Instead, in line 14, we re-privatize the private model with much
lower noise of variance $s\sigma_{2}^{2}$. We have found that this re-
privatization step can significantly enhance the performance of our proposed
algorithm, ensuring its sparse basis recovery success. Please see the
numerical results in Section 7 for details. Intuitively, although we re-
privatize the same features/artifacts twice, since, we add lower noise in the
second time, the noise in $\tilde{\alpha}$ estimate (refer line (15) in
Algorithm 3) is of order only depending on $s$. In contrast, if we do not
utilize this modification, then the order of noise in both $\bm{\beta}$ and
$\tilde{\bm{\alpha}}$ will be $\mathcal{O}(\sqrt{p})$. Since these estimates
will be used in later steps to identify the next basis, without re-
privatization the higher noise will disrupt the correlation computation,
leading to stricter requirements on the sample size.

###  3.1 The SPriFed-OMP-GRAD Algorithm

We now present the second version of private OMP with our second enhancement
in Algorithm 4 that further enhances the performance of Algorithm 3.

_Enhancement 2: Re-visiting SPriFed-OMP from the gradient perspective_

In line 15 of Algorithm 3, we compute the private correlation (essentially the
gradient) at the $l^{th}$ step given by,
$\bm{\gamma_{0}^{-}-\beta\beta_{S_{\mathcal{A}}}\gamma_{S_{\mathcal{A}}}=(\gamma_{0}^{-})^{t}+\eta_{\gamma_{0}}-(\beta^{t}+\eta_{\beta})\tilde{\alpha}}$
where $\bm{(\gamma_{0}^{-})^{t},\beta^{t}}$ represent the true (non-noisy
correlations), $\bm{\eta_{\gamma_{0}},\eta_{\beta}}$ represent the
corresponding noise values and
$\bm{\tilde{\alpha}(=\beta_{S_{\mathcal{A}}}\gamma_{S_{\mathcal{A}}})}$
represents the privatized linear model. We notice that since we are already
privatizing the correlations required in the computation of
$\bm{\tilde{\alpha}}$, we are naturally inclined also to privatize
$\bm{(\gamma_{0}^{-})^{t}}$ and $\bm{\beta^{t}}$ similarly. Thus, in our first
approach of Algorithm 3, all clients collaboratively compute correlations, and
then the server computes the final private model and the private gradient.
However, our second approach improves this correlation computation by first
computing the $L_{2}$ sensitivity (Definition A.2) of the gradient calculated
using the private model (or mathematically,
$\bm{(\gamma_{0}^{-})^{t}-\beta^{t}\tilde{\alpha}}$) and then adding the
corresponding DP noise to the entire correlation rather than separately
privatizing each correlation. In this second approach, given the model, each
client computes the gradient on its own device. Collaboratively, all clients
and the server then compute the aggregated gradient. However, the server still
handles the model computation with the noisy correlation method. That way, the
server can share the private model with individual clients without sacrificing
privacy. Later in the next section and the experimental section, we will also
discuss a clipping-based version of the new algorithm, SPriFed-OMP-GRAD
(Algorithm 4) that leverages clipping to provide better empirical performance.

We now present the second major enhancement in our proposed method in
Algorithm 4. We note that the correlation computed in line 15 of Algorithm 3
is essentially the private gradient computed at the current model value. In
Algorithm 3, we first separately privatize correlation values and then combine
them to obtain the private gradient. On the other hand, in Algorithm 4, we let
the server pass the private model to the clients, and let the clients directly
compute the residues and the local gradients (_i.e.,_ the gradient computed on
line 4 in Algorithm 4). Then, we aggregate the local gradients in a private
manner. We expect that Algorithm 4 will see a significant performance benefit
due to the effect of clipping. Note that clipping can control the sensitivity
of individual data items, and thus has a directly impact on the magnitude of
DP noise needed. Since, gradients typically reduce in magnitude as training
proceeds, we can expect to clip the gradients more aggressively without
affecting performance. In contrast, the correlations in Algorithm 3 will stay
at large values during the entire training, which is not amenable to clipping.
We report experimental effects of clipping on gradients and correlation in
section 7 to visualize their varying impact.

Algorithm 1 NOISY-SMPC

1:Input Parameters: Vectors to be multiplied: $\bm{X_{k}},\bm{X_{j}}\text{ (or
}\bm{y})\in\mathbb{R}^{n}$, Number of clients: $n$, Noise Variance:
$\sigma_{0}^{2}$

2:Output:
$\bm{X_{k}^{T}X_{j}+\eta_{*}};\bm{\eta_{*}}\sim\mathcal{N}(0,\sigma^{2}_{0})$

3:procedure NOISY-SMPC($\bm{X_{k}},\bm{X_{j}},n,\sigma_{0}$)

4: Server broadcasts noise variance $\sigma_{0}^{2}$ to all clients

5: for client $i\in[n]$ do

6: $\eta_{i}\sim\mathcal{N}(0,\frac{\sigma_{0}^{2}}{n}$)

7: Compute $q_{i}\leftarrow\bm{X}_{k,i}\bm{X}_{j,i}+\eta_{i}$

8: end for

9: // Vanilla SMPC

10: Compute $q_{total}\leftarrow\sum_{i=1}^{n}q_{i}$ using SMPC (Bonawitz et
al., 2017)

11: // The resulting private sum satisfies the following $\rightarrow
q_{total}=\sum_{i=1}^{n}\bm{X}_{k,i},\bm{X}_{j,i}+\eta_{i}=\bm{X_{k}^{T}}\bm{X_{j}}+\sum_{i=1}^{n}\eta_{i}=\bm{X_{k}^{T}}\bm{X_{j}}+\eta_{*}$;
$\eta_{*}\sim\mathcal{N}(0,n\cdot\frac{\sigma_{0}^{2}}{n})(=\mathcal{N}(0,\sigma_{0}^{2}))$

12: Server obtains and returns $q_{total}$ to all clients

13:end procedure

Algorithm 2 PRIVATE-OLS

1:Input Parameters: Data:
$\bm{X}=\\{\bm{X_{1},X_{2},...,X_{p}}\\}\in\mathbb{R}^{n\times p}$, where
$\bm{X_{i}}$ is the $i^{th}$ feature. Response Variable: $y\in\mathbb{R}^{n}$,
Feature Set: $S_{\mathcal{A}}$, New Feature to be privatized: $l^{*}$, Privacy
Parameter: $\sigma_{2}$, True Model Support Cardinality: $s$.

2:Output: Private Model:
$\tilde{\alpha}=(X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}+\eta_{\beta_{S_{\mathcal{A}}}})^{-1}(X_{S_{\mathcal{A}}}^{T}y+\eta_{2});\eta_{\beta_{S_{\mathcal{A}}}}$
and $\eta_{\gamma_{S_{\mathcal{A}}}}$ both have iid elements with distribution
$\mathcal{N}(0,\sigma^{2}_{2}s)$.

3:procedure PRIVATE-OLS($\bm{X},\bm{y},S_{\mathcal{A}},l^{*},\sigma_{2},s$)

4: $\bm{\gamma_{S_{\mathcal{A}}}}[l]\leftarrow\textbf{NOISY-
SMPC}(\bm{X_{l^{*}}},y,n,\sigma_{2}\sqrt{s})$ /* Note, that total private
correlation over all previous rounds is
$\bm{\gamma_{S_{\mathcal{A}}}=X_{S_{\mathcal{A}}}^{T}y+\eta_{\gamma_{S_{\mathcal{A}}}}};\bm{\eta_{\gamma_{S_{\mathcal{A}}}}}\in\mathbb{R}^{l\times
1}$ has iid elements from $\mathcal{N}(0,\sigma_{2}^{2}s)$ */

5: for $k=0:l$ do

6: $\bm{\beta_{S_{\mathcal{A}}}}[l,k]\leftarrow\textbf{NOISY-
SMPC}(\bm{X_{l^{*}},X_{k}},n,\sigma_{2}\sqrt{s})$

7: end for /* Note that, as $\bm{\beta_{S_{\mathcal{A}}}}$ is symmetric, we
have $\bm{\beta_{S_{\mathcal{A}}}}[k,l]=\bm{\beta_{S_{\mathcal{A}}}}[l,k]$ and
thus total private covariance is
$\bm{\beta_{S_{\mathcal{A}}}\triangleq\widehat{X^{T}{S_{\mathcal{A}}}X_{S_{\mathcal{A}}}}+\eta_{\beta_{S_{\mathcal{A}}}}};\bm{\eta_{\beta_{S_{\mathcal{A}}}}}\in\mathbb{R}^{s\times
s}$ has iid elements from $\mathcal{N}(0,\sigma_{2}^{2}s)$ */

8: Return $\bm{(\beta_{S_{\mathcal{A}}})^{-1}\gamma_{S_{\mathcal{A}}}}$

9:end procedure

Algorithm 3 SPriFed-OMP: Private Orthogonal Matching Pursuit

1:Input Parameters: Data:
$\bm{X}=\\{\bm{X_{1},X_{2},...,X_{p}}\\}\in\mathbb{R}^{n\times p}$, where
$\bm{X_{i}}$ is the $i^{th}$ feature. Response Variable: $y\in\mathbb{R}^{n}$,
Feature Set: $\Omega\triangleq\\{1,2,...,p\\}$, True Model Support
Cardinality: $s$, GDP Privacy Parameters: $\mu_{p},\mu_{s}(\mu_{p}>\mu_{s})$,
Noise standard deviation:
$\sigma_{1}=\frac{1}{\mu_{p}},\sigma_{2}=\frac{1}{\mu_{s}}$.

2:Output: Predicted Support $S_{\mathcal{A}}$, Predicted sparse model
$\hat{\alpha}$ over predicted support $S_{\mathcal{A}}$

3:Initialize: $S_{\mathcal{A}}=\emptyset,\bm{X_{S_{\mathcal{A}}}}=\bar{0},l=0$

4:for $k=0:p-1$ do /* Server and Clients privately compute data-response
correlation */

5: $(\bm{\gamma_{0}^{-}})[k]\leftarrow\textbf{NOISY-
SMPC}(\bm{X_{k}},\bm{y},n,\sigma_{1}\sqrt{p})$

6:end for /* Note that the total private correlation overall $k$ is
$\bm{\gamma_{0}^{-}\triangleq\widehat{X^{T}y}=X^{T}y+\eta_{\gamma}};\bm{\eta_{\gamma}}\in\mathbb{R}^{p\times
1}$ has iid elements from $\mathcal{N}(0,\sigma_{1}^{2}p)$ */

7:while $|S_{\mathcal{A}}|\leq s$ do

8: $l^{*}\leftarrow\operatorname*{arg\,max}_{j\in\Omega/S_{\mathcal{A}}}\text{
}|(\gamma_{l}^{-})_{j}|$ /* Server privately extracts new highest-correlated
feature */

9: $S_{\mathcal{A}}\leftarrow S_{\mathcal{A}}\cup{l^{*}}$

10: for $k=0:p-1$ do /* Server and Clients privately compute data covariance
for newly extracted feature */

11: $\bm{\beta}[k,l]\leftarrow\textbf{NOISY-
SMPC}(\bm{X_{k},X_{l^{*}}},n,\sigma_{1}\sqrt{p})$

12: end for /* Note that, the total private covariance over all previous
rounds is  $\bm{\beta=\widehat{X^{T}X_{S_{\mathcal{A}}}}\triangleq
X^{T}X_{S_{\mathcal{A}}}+\eta_{\beta}};\bm{\eta_{\beta}}\in\mathbb{R}^{p\times(l+1)}$
has iid elements from $\mathcal{N}(0,\sigma_{1}^{2}p)$ */

13: // Server privately updates residual with newly chosen feature

14: $\tilde{\alpha}\leftarrow$ PRIVATE-
OLS($\bm{X,y},S_{\mathcal{A}},l^{*},\sigma_{2},s$)

15: $\bm{\gamma_{l}^{-}\leftarrow\gamma_{0}^{-}-\beta\tilde{\alpha}}$

16: $l=l+1$

17:end while

18:Return $S_{\mathcal{A}}$,
$\bm{\hat{\alpha}\leftarrow(\beta_{S_{\mathcal{A}}})^{-1}\gamma_{S_{\mathcal{A}}}}$
/* Server privately publishes sparse basis and model overall extracted
features */

Algorithm 4 SPriFed-OMP-GRAD: Private Orthogonal Matching Pursuit

1:Input Parameters: Data:
$\bm{X}=\\{\bm{X_{1},X_{2},...,X_{p}}\\}\in\mathbb{R}^{n\times p}$, where
$\bm{X_{i}}$ is the $i^{th}$ feature. Response Variable: $y\in\mathbb{R}^{n}$,
Feature Set: $\Omega\triangleq\\{1,2,...,p\\}$, True Model Support
Cardinality: $s$, GDP Privacy Parameters: $\mu_{p},\mu_{s}(\mu_{p}>\mu_{s})$,
Noise standard deviation:
$\sigma_{1}=\frac{1}{\mu_{p}},\sigma_{2}=\frac{1}{\mu_{s}}$,  Privacy
Constants:
$X_{M}=y_{M}=1,B_{C}=1+2X_{M}\kappa_{\varepsilon}+2\sqrt{s}X_{M}^{2}(\frac{\lVert(\alpha_{*})_{S^{c}_{\mathcal{A}}}\rVert_{\infty}}{1-\zeta_{s+1}}+\frac{\sqrt{1+\zeta_{s+1}}\kappa_{\varepsilon}}{1-\zeta_{s+1}})$
(refer Lemma 5 and Section 7 for further details about setting $B_{C}$ in
practice).

2:while $\lvert S_{\mathcal{A}}\rvert\leq s$ do

3: for $j=0:p$ do

4: $\tilde{\gamma_{l}}[j]\leftarrow\textbf{NOISY-
SMPC}(\bm{X_{j}},\bm{y-X_{S_{\mathcal{A}}}\tilde{\alpha}},n,\sigma_{1}B_{C}\sqrt{p})$

5: end for

6:
$l^{*}\leftarrow\operatorname*{arg\,max}_{j\in\Omega/S_{\mathcal{A}}}\lvert(\tilde{\gamma_{l}})_{j}\rvert$

7: $S_{\mathcal{A}}\leftarrow S_{\mathcal{A}}\cup l^{*}$

8: // Server privately updates model with newly chosen feature

9: $\tilde{\alpha}\leftarrow$ PRIVATE-
OLS($\bm{X,y},S_{\mathcal{A}},l^{*},\sigma_{2},s$)

10: $l=l+1$

11:end while

##  4 Privacy Analysis

This section includes the details for the privacy cost incurred by running
Algorithms 3 and 4.

Since we need to account for the overall DP guarantee by composing multiple DP
mechanisms, in this paper, we will use Gaussian Differential Privacy (GDP)
(Dong et al., 2019) (a variant of DP). We choose GDP for its superior
composition performance compared to other methods, such as advanced
composition (Kairouz et al., 2015), Renyi DP (Mironov, 2017), and concentrated
DP (Bun & Steinke, 2016). Further, it can be easily converted to DP. Before
stating the Gaussian mechanism, we first formally define the
$L_{2}$-sensitivity of a function.

######  Definition 4.1 ($L_{2}$ sensitivity of a function).

Given a deterministic function $f$, consider all possible pairs of neighboring
datasets $\bm{X}$ and $\bm{X^{{}^{\prime}}}$ in the domain of $f$ that differ
on a single row. Then, we denote $l_{2}$-sensitivity of $f$ as

$\displaystyle\Delta_{2}f\triangleq\max_{\bm{X},\bm{X^{{}^{\prime}}}\in
dom(f)}||f(\bm{X})-f(\bm{X^{{}^{\prime}}})||_{2}$

Thus, $\Delta_{2}f$ records the maximum possible change in $f$ due to a change
in a single row. Often, each row belongs to a different client, and thus,
$\Delta_{2}f$ will record the maximum change due to a change in a client or
due to the inclusion or removal of a client from the dataset. Below, we state
the Gaussian mechanism from Dong et al. (2019) required to achieve the GDP
guarantee.

######  Lemma 1.

[GDP mechanism (Theorem 2.7 from Dong et al. (2019))] Let
$f:\mathbb{R}^{n\times p}\rightarrow\mathbb{R}^{p}$ be some function computed
over the dataset $\bm{X}\in\mathbb{R}^{n\times p}$. Then, the randomized
Gaussian mechanism $\mathcal{M}(\bm{X})=f(\bm{X})+\eta$ is $\mu-GDP$ where
$\eta\sim\mathcal{N}(0,\frac{(\Delta_{2}f)^{2}}{\mu^{2}}\mathbb{I}_{p})$ and
$\Delta_{2}f$ is the $L_{2}$\- sensitivity of $f$ (Definition A.2).

Note that GDP guarantees can be easily and optimally composed over multiple
iterations via Corollary 3.3 in Dong et al. (2019). Furthermore, we can
convert between a $\mu-GDP$ guarantee and the $(\varepsilon,\delta)$-DP
guarantee losslessly via Corollary 2.13 in Dong et al. (2019). However, we
only need to convert $\mu$-GDP to DP for our purpose. Below, we present both
of these lemmas.

######  Lemma 2.

[Conversion from $\mu$-GDP to $(\varepsilon,\delta)$-DP] (Corollary 2.13 in
Dong et al. (2019)) A mechanism is $(\varepsilon,\delta(\varepsilon))$-DP for
all $\varepsilon\geq 0$ and

$\displaystyle\delta(\varepsilon)=\Phi(-\frac{\varepsilon}{\mu}+\frac{\mu}{2})-e^{\varepsilon}\Phi(-\frac{\varepsilon}{\mu}-\frac{\mu}{2}).$

if it is $\mu$-GDP. Here, $\Phi$ denotes the standard normal CDF.

Further, multiple mechanisms with $\mu$-GDP into an overall GDP mechanism.

######  Lemma 3.

Composition of $\mu$-GDP (Corollary 3.3 in Dong et al. (2019)) The composition
$k$-GDP mechanisms, each of which is $\mu_{i}$-GDP, $i=1,...,k$ is
$\sqrt{\sum_{i=1}^{k}{\mu_{i}^{2}}}$-GDP.

In summary, leveraging GDP allows us to obtain better composition results,
control the privacy leakage by a single parameter, and losslessly convert
between GDP and $(\varepsilon,\delta)$-DP. For further details regarding GDP,
we refer the reader to Dong et al. (2019).

For our proposed Algorithm 3, we notice that we only need to differentially
privatize the terms $\bm{X_{j}^{T}y}$ and $\bm{X^{T}X_{j}}$ and their
combinations. Thus, we use Lemma 4 to identify their $L_{2}$-sensitivity
(Definition A.2).

######  Lemma 4.

Under Assumption 2, the $L_{2}$-sensitivity (Definition A.2) of the terms
$\bm{X_{j}^{T}y}$ and $\bm{X^{T}X_{j}}$ for any $j\in S_{\mathcal{A}}$ as
defined in Algorithm 3 are given by
$2\sqrt{p}X_{M}y_{M}=\mathcal{O}(\sqrt{p})$ and
$2\sqrt{p}X_{M}^{2}=\mathcal{O}(\sqrt{p})$ respectively. Here,
$\bm{X}\in\mathbb{R}^{n\times p}$ and $\bm{y}\in\mathbb{R}^{n}$.

###### Proof.

The proof is available in the Appendix section C. ∎

######  Lemma 5.

Consider a design matrix $X\in\mathbb{R}^{n\times p}$ satisfying RIP of order
$s+1$ and RIC $\zeta_{s+1}$. Now given the correlation/gradient $C_{j}$ (as
computed in Line 7 of algorithm 4) for any $j\in[p]$, we show that the upper
bound on the $L_{2}$ sensitivity of the expression is upper bounded by the
following values,

$\displaystyle\Delta_{2}(C_{j})\leq
1+2X_{M}\kappa_{\varepsilon}+2\sqrt{s}X_{M}^{2}(\frac{\lVert(\alpha_{*})_{S^{c}_{\mathcal{A}}}\rVert_{\infty}}{1-\zeta_{s+1}}+\frac{\sqrt{1+\zeta_{s+1}}\kappa_{\varepsilon}}{1-\zeta_{s+1}})\triangleq
B_{C}$

when Assumptions 1 and 2 hold. Here, $\alpha_{*}$ is the underlying ground-
truth model, $\kappa_{\varepsilon}=\sqrt{2\log(2n/p_{b})}\sigma_{\varepsilon}$
and $\sigma_{\varepsilon}$ is the additive system error’s standard deviation.
We also need to assume that $\zeta_{s+1}<\frac{1}{\sqrt{s}}$ and
$n>\frac{6s^{2}C_{1}X_{M}^{5}y_{M}\sqrt{2\log(2s/p_{b})}}{\mu_{s}^{2}(1-\zeta_{s+1})^{2}}$.
The result holds with probability $1-3p_{b}$ where $p_{b}$ is a small positive
probability value.

###### Proof.

The proof is provided in the Appendix Section F. ∎

######  Theorem 6.

Algorithms 3 and 4 satisfies $\mu$-GDP with
$\mu=\sqrt{s\cdot\mu_{p}^{2}+2s\cdot\mu_{s}^{2}}$, where $\mu_{p},\mu_{s}$ are
privacy constants as inputs to Algorithm 3.

###### Proof.

The proof follows simply by Lemma 3. In our proposed algorithm, we use the
$\mu_{p}$-GDP (lines 5 and 11 in Algorithm 3 and line 4 in Algorithm 4)
mechanism $s$ times and the $\mu_{s}$-GDP mechanism $2s$ times (lines 4 and 6
in algorithm 2). Thus, the result follows. Finally, using Lemma 2, we can
obtain the corresponding DP guarantee of Algorithm 3. ∎

##  5 Accuracy of Private Orthogonal Matching Pursuit

This section will present our main result on the accuracy of sparse basis
recovery. Theorem 7 states that the proposed algorithm SPriFed-OMP (Algorithm
3) terminates in finite steps and can recover the true sparse basis and the
corresponding model w.h.p. as $n\rightarrow\infty$ and when
$n=\mathcal{O}(\sqrt{p})$.

Before we state our main results, we define the relevant constants. Let
$\zeta_{s+1}$ denote the RIC of $\bm{X}$ of $(s+1)$-order. Let $\nu$ be a
small positive constant, close to $1$, that satisfies the inequality
$\frac{1}{1-\zeta_{s+1}}\leq(1+\nu\zeta_{s+1})$. Next, we define
$\kappa_{\epsilon}=\sigma_{\epsilon}\sqrt{2\log(\frac{2n}{p_{b}})}$,
$\kappa_{1}=\kappa_{32}(1+\nu\zeta_{s+1})^{3/2}$ and
$\kappa_{2}=\kappa_{1}(1+\nu\zeta_{s+1})^{5/2}(1+(\nu+1)\zeta_{s+1})$ where
$\kappa_{32}=\kappa_{M_{2}}\sqrt{\log(p/p_{b})}/\mu_{p}$,
$\kappa_{1}=\kappa_{M_{2}}\sqrt{\log(s/p_{b})}/\mu_{s}$. Here, $p_{b}$ is a
small positive probability, eventually affecting the with-high-probability
statement of Theorem 7. Further define the constant $\kappa_{M_{2}}$ as

$\displaystyle\kappa_{M_{2}}$
$\displaystyle=1+\frac{C_{1}s^{3/2}}{n\mu_{s}(1-\zeta_{s+1})},$

where the constant $C_{1}$ (which depends on $p_{b}$) is chosen such that the
following probability statement holds: For a square matrix $N$ with dimensions
$s\times s$, where the elements follow an i.i.d. distribution of
$\mathcal{N}(0,1)$, its largest eigenvalue is no greater than $C_{1}\sqrt{s}$
with probability $1-p_{b}$, _i.e.,_

$\displaystyle\Pr\Big{(}||N||_{2}>C_{1}\sqrt{s}\Big{)}\leq p_{b}$

By these choices (which will be used in our proof in the technical report),
each constant bound the corresponding random variable with probability
$1-p_{b}$. Finally, we denote the true basis of the underlying model by
$S_{*}$.

######  Theorem 7.

[Sparse recovery of SPriFed-OMP] Under the provided system model (Section 2)
and Assumptions 1 and 2, suppose that the following conditions holds:

  1. 1.

Sample size:
$n\geq\max\\{4\sqrt{p}s\kappa_{1},16s^{5/2}\kappa_{2},\frac{4\sqrt{log(p-s)}\sqrt{p}}{\mu_{p}\lVert\alpha_{S_{*}}\rVert_{2}}\\}=\mathcal{O}(\max\\{s\sqrt{p\log
p},s^{5/2}\\})$

  2. 2.

RIC of $\bm{X}$ of order-$s+1$: $\zeta_{s+1}\leq\frac{1}{4(1+\sqrt{s})}$

  3. 3.

$\kappa_{\epsilon}$ : $\kappa_{\epsilon}\leq\frac{min_{j\in
S_{*}}|\bm{\alpha_{j}}|}{16(\sqrt{s}+1)(1+\nu\zeta_{s+1})}$

Then, the SPriFed-OMP algorithm will correctly recover the true basis with
probability $1-40(s+1)\cdot p_{b}$.

_Remark (Intuition behind the assumptions for successful recovery in Theorem
7)_: If we ignore the $\log{p}$ term, condition 1 only requires
$n=\mathcal{O}(\sqrt{p})$ for sparse recovery. To the best of our knowledge,
this is the first result in the literature that both attains DP in an FL
setting and attains sparse recovery when $n$ is much smaller than $p$. We note
that condition 1 is stricter than typical OMP results in the non-private
settings, where $n=\mathcal{O}(\log{p})$ is sufficient. We believe the main
reason for this difference is that we must simultaneously ensure DP in the FL
setting. As we discussed earlier when we present Algorithm 3, even in the
first step for computing the data-response correlation,
$\mathcal{O}(\sqrt{p})$ noise seems unavoidable (remark under Section 3).
Thus, it seems difficult to suppress the noise when $n$ is smaller than
$\sqrt{p}$.

Condition 2 is related to the strength of the RIP assumption. The RIC is an
important attribute of the RIP. The smaller the RIC, the closer $\bm{X}$ is to
a unitary matrix, and thus, the corresponding condition 2 becomes more
demanding. Our condition 2 provides a simple upper bound for the RIC. As
expected, we observe that, with the increasing value of $s$, our RIC bound
reduces, implying that sparse recovery gets more challenging as the model
cardinality increases. This intuition is in line with that in the non-private
OMP literature.

Condition 3 on the additive error in the training data is also intuitive. If
the minimum true-model column norms are large enough (_i.e.,_ the signals due
to true features are sufficiently strong), then we expect that the impact of
the additive error in the system will be minuscule. This is reflected by the
term $\min_{j\in S_{*}}|\bm{\alpha_{j}}|$ in the numerator. We also expect
that model recovery gets more complex with rising model cardinality (_i.e.,_
rising $s$). Thus, the impact of the error in the training data will rise as
reflected by the $\sqrt{s}+1$ term in the denominator. Similar to Theorem 8,
we can state the recovery guarantee of SPriFed-OMP-GRAD below.

######  Theorem 8.

[Sparse Recovery of SPriFed-OMP-Grad] Under the provided system model (Section
2) and Assumptions 1 and 2 suppose that the following conditions hold,

  1. 1.

Sample Size:
$n\geq\max\Bigg{\\{}\frac{4(\sqrt{s}+1)(1+\sqrt{s}(\kappa_{\alpha}+1))\sqrt{p}\kappa_{p}^{{}^{\prime}}}{\mu_{p}\lVert\alpha_{S_{\mathcal{A}}^{c}}\rVert_{2}},\frac{4\sqrt{\log(p-s)}\sqrt{p}}{\mu_{p}\lVert\alpha_{S_{*}}\rVert_{2}}\Bigg{\\}}=\mathcal{O}(s\sqrt{p\log
p})$

  2. 2.

RIC of $X$ of order $s+1$: $\zeta_{s+1}\leq\frac{1}{\sqrt{s}+1}$

  3. 3.

$\kappa_{\varepsilon}<\frac{\min_{j\in
S_{*}}\lvert\alpha_{j}\rvert}{6(1+\nu\zeta_{s+1})^{3/2}\Big{(}\sqrt{s}+1+\frac{C_{1}\sqrt{s}\kappa_{s}^{{}^{\prime}}\sqrt{1+\zeta_{s+1}}}{n\mu_{s}}\Big{)}}$

Then the SPriFed-OMP-GRAD algorithm will correctly recover the true basis with
probability $1-8(s+1)p_{b}$. Note that here,
$\kappa_{\alpha}=\frac{\sqrt{1+\zeta_{s+1}}(\sqrt{1+\zeta_{s+1}}\lVert\alpha_{*}\rVert+\kappa_{\varepsilon})}{1-\zeta_{s+1}},\kappa_{p}^{{}^{\prime}}=\sqrt{2\log(\frac{2p}{p_{b}})},\kappa_{s}^{{}^{\prime}}=\sqrt{2\log(\frac{2s}{p_{b}})}$
and
$B_{C}=1+2X_{M}\kappa_{\varepsilon}+2\sqrt{s}X_{M}^{2}(\frac{\lVert(\alpha_{*})_{S^{c}_{\mathcal{A}}}\rVert_{\infty}}{1-\zeta_{s+1}}+\frac{\sqrt{1+\zeta_{s+1}}\kappa_{\varepsilon}}{1-\zeta_{s+1}})$
and $p_{b}$ is a small positive probability value similar to Theorem 7.

Remark (Benefits of Algorithm 4 over Algorithm 3: Readers can see that the
conditions of Theorem 8 are on the same order as that of Theorem 7. However,
this is because we assume the DP noise magnitudes are the same in both cases.
In practice and as we discussed earlier, due to the decreasing magnitude of
the gradients, Algorithm 4 can benefit from more aggressive clipping and lower
DP noise magnitude. We will present numerical results to demonstrate the
benefits of these enhancement in Section 7.

Based on our main Theorems 7 and 8, next, we provide results on the parameter
estimation error and empirical risk. Both of these results are optimal as they
are in the same order as several previous results for both the parameter
estimation error (Meinshausen & Yu, 2009; Wasserman & Roeder, 2009) and the
empirical risk analysis (Kifer et al., 2012; Bassily et al., 2014) even under
the assumption that they have already known the correct basis before-hand. In
particular, these error bounds depend on $s$ but are independent of the
dimension $p$ since we have already identified the correct sparse basis as in
Theorems 7 and 8.

######  Theorem 9.

[Estimation Error of SPriFed-OMP and SPriFed-OMP-Grad] Under the same system
model and assumptions as Theorem 7, with a high probability of $1-p_{b}$, the
estimation error satisfies,

$\displaystyle\Delta\bm{\alpha}$
$\displaystyle\triangleq||\bm{\hat{\alpha}-\alpha_{S_{*}}}||_{2}$
$\displaystyle\leq\frac{(1-\zeta_{s+1})\sqrt{2s\log(2s/p_{b}))}\sigma_{\epsilon}X_{M}}{\sqrt{n}}+\frac{s\sqrt{2\log(\frac{2s}{p_{b}})}}{n\mu_{s}(1-\zeta_{s+1})}$
$\displaystyle+\frac{s^{3/2}\kappa_{M}\sqrt{2\log(\frac{2s^{2}}{p_{b}})}}{\mu_{s}n(1-\zeta_{s+1})^{2}}\Big{(}(1+\zeta_{s+1})||\alpha_{S_{*}}||_{2}$
$\displaystyle+\sqrt{1+\zeta_{s+1}}\kappa_{\epsilon}+\frac{\sqrt{s}\sqrt{2\log(\frac{2s}{p_{b}})}}{\mu_{s}n}\Big{)}$
$\displaystyle=\mathcal{O}\Big{(}\sqrt{\frac{s\log(s)}{n}}\Big{)}$

where $\bm{\hat{\alpha}}$ is the predicted model from the SPriFed-OMP
algorithm, $\bm{\alpha_{S_{*}}}$ is the true ground-truth model, $\zeta_{s+1}$
is the RIC constant for the matrix $\bm{X}$ and $\mu_{s}$ is as defined in
Algorithm 3. Our result holds with a high probability of $1-4p_{b}$.

Proof: The proof is omitted due to page limits.

######  Theorem 10.

[Risk Analysis of SPriFed-OMP] Under the same system model and assumptions as
Theorem 7, with probability $1-4p_{b}$ the risk $\Delta R$ from algorithm
SPriFed-OMP is bounded by,

$\displaystyle\Delta R$ $\displaystyle\triangleq
R(\bm{\hat{\alpha}};\bm{X},\bm{y},n)$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\Big{(}(\bm{x_{i}}\bm{\hat{\alpha}}-y_{i})^{2}-(\bm{x_{i}}\bm{\alpha_{S_{*}}}-y_{i})^{2}\Big{)}$
$\displaystyle\leq
2\Big{(}\frac{(1+\zeta_{s+1})s\sqrt{2\log(\frac{2s}{p_{b}})}}{\mu_{s}n(1-\zeta_{s+1})}\Big{)}^{2}$
$\displaystyle+2\Big{(}\frac{(1+\zeta_{s+1})\sqrt{2ns\log(2s/p_{b}))}\sigma_{\epsilon}X_{M}}{n(1-\zeta_{s+1})}\Big{)}^{2}$
$\displaystyle+2\Big{(}\frac{s^{3/2}\kappa_{M}\sqrt{2\log(\frac{2s^{2}}{p_{b}})}}{\mu_{s}n(1-\zeta_{s+1})^{2}}\Big{(}(1+\zeta_{s+1})||\alpha_{S_{*}}||_{2}$
$\displaystyle+\sqrt{1+\zeta_{s+1}}\kappa_{\epsilon}+\frac{\sqrt{s}\sqrt{2\log(\frac{2s}{p_{b}})}}{\mu_{s}n}\Big{)}\Big{)}^{2}$
$\displaystyle=\mathcal{O}\Big{(}\frac{2s\log(s))}{n}\Big{)}$

where $\bm{\hat{\alpha}}$ is the predicted model from the SPriFed-OMP
algorithm, $\bm{\alpha_{S_{*}}}$ is the true ground-truth model, $\zeta_{s+1}$
is the RIC constant for the matrix $\bm{X}$ and $\mu_{s}$ is as defined is
defined in Algorithm 3. $\bm{X,y}$ are the data matrix and the label vector
while $n$ is the number of clients.

Proof: The proof is omitted due to page limits.

##  6 Proof Sketch of Theorems 7 and 8

This section identifies the key steps required to prove Theorems 7 and 8. From
Algorithms 3 and 4, we notice that in the very first round, the server picks
the feature with a maximum absolute data-response correlation over the entire
feature set. To ensure that we only pick features from the true feature set,
we thus require expression 1 below to hold.

$\displaystyle\min_{j\in
S_{*}}|\bm{X_{j}^{T}y+\eta_{\gamma_{j}}}|\geq\max_{j\notin
S_{*}}|\bm{X_{j}^{T}y+\eta_{\gamma_{j}}}|.$ (1)

In other words, for correct basis recovery to occur in the first round, the
minimum absolute feature correlation over the true set should exceed the
maximum absolute feature correlation over all the features not in this true
set. The strategy for proving 1 is to bound the randomness-inducing terms
(_i.e.,_ noise terms) with high probability. This high-probability event will
also dictate the relationship between the sample size, the restricted isometry
constant, and the rest of the system constants. Using the restricted isometry
property over the dataset contributed by all clients, we can then show that
the above inequality holds. The remaining rounds are handled similarly over
the updated residuals. The complete proof is available in the Appendix
sections E and F.

Remark: Here, we re-iterate that both algorithms’ analyses have subtle but
critical differences. In Algorithm 3, the noise is added to column
correlations (that are collaboratively computed by the clients), and then the
private model and the aggregate correlation between columns and the residual
(_i.e.,_ the gradient) are both computed by the server. In Algorithm 4, the
clients first compute their individual gradients and then collaboratively find
the average private gradient. The server computes the private model similarly
to in Algorithm 3 but then shares this model with the individual clients.
Here, the clients then recompute their gradients for the next step. Algorithm
3 adds the noise to the correlations of the gradient, which are then
multiplied. On the other hand, in Algorithm 4, the bulk of the noise value is
added outside of the correlations, directly to the gradient.

Although our analysis idea is similar to several OMP proofs in the literature
(such as Chen & Caramanis (2012)), we have to redo the analysis because the
way with which we add noise to the system (to ensure DP) is completely
different. Notably, Chen & Caramanis (2012) adds noise to the measurement
matrix $\bm{X}$, which has several issues. First, it does not even protect the
privacy of the response $y$. Second, the noise variance is of the order-$np$,
which is too high and likely leads to poor estimation accuracy. Indeed, the
estimation error bound (Corollary 8) in Chen & Caramanis (2012) is given by
$||\bm{\alpha_{*}}-\bm{\hat{\alpha}}||_{2}=\mathcal{O}(\frac{(\sigma_{w}+\sigma_{w}^{2})||\bm{\alpha_{*}}||_{2}\sqrt{s\log
p}}{n})$, where $\sigma_{w}$ is the standard deviation of the noise matrix
$\bm{W^{T}W}$ such that $\bm{W}$ is the additive noise added to the original
matrix $\bm{X}$. From their analysis, the term
$\sigma_{w}+\sigma_{w}^{2}=\mathcal{O}(p^{2})$ when $\sigma_{w}$ is set to the
DP-compliant noise value. The above error bound, therefore, will be huge when
$p\gg n$, negating its usefulness in the high-dimensional regime. The
experimental results in Chen & Caramanis (2012) assume small values of
$\sigma_{w}$, not at the level of $\mathcal{O}(p)$ needed for achieving DP.
Third, Chen & Caramanis (2012) does not provide a condition on the RIC. They
assume that the matrix elements are subgaussian i.i.d.; thus, their model
implicitly determines the RIC constant in their results. In contrast, we allow
an arbitrary design matrix $\bm{X}$ (irrespective of the matrix’s data
distribution) that satisfies an RIP condition. Thus, it is important to
quantify the condition of RIC to guarantee sparse recovery.

##  7 Empirical Results

In this section, we compare our proposed algorithm SPriFed-OMP with a version
of the DP-SGD algorithm (Abadi et al., 2016) that uses $L_{1}$ regularization
(Zhao & Yu, 2006). The first set of experiments is based on a synthetic data
set so that we can freely vary the sample size $n$, the model dimension $p$,
and the model sparsity $s$. Given $p$ and $n$, we first generate $\bm{X}$ with
_i.i.d._ elements, each of which is from $\mathcal{N}(0,1)$. Given the model
sparsity $s$, we then generate the non-zero coefficients of the model
parameters $\bm{\alpha_{*}}$, each of which is from the distribution
$\mathcal{N}(2,1)$. We then add additive error with mean zero and standard
deviation $\sigma_{\epsilon}=0.001$ to generate the response $y$. We then clip
(to $1$ or $-1$) every element in $\bm{X}$ and $\bm{y}$ with magnitude greater
than $1$ and then re-scale the whole matrix/vector appropriately to maintain
their unit variance as required by the RIP condition 1.

###  7.1 Warm-up Experiments: Intuition behind our proposed enhancements

Before moving on to the main experimental results, we revisit the enhancements
noted in Section 3. Note that both experiments below use the synthetic data
setup above.

Importance of Enhancement 1: Adding lower noise to selected features

Here, we test the sparse basis recovery capabilities of Algorithm 3 on a
synthetic toy dataset. We compare Algorithm 3 with another modified version,
where Enhancement 1 is removed, i.e., we directly use the artifacts privatized
by the higher order noise in lines 5 and 11 to estimate the private model.
Thus, the private model will now have noise of order $\mathcal{O}(\sqrt{p})$
instead of $\mathcal{O}(\sqrt{s})$. In Table 1, we observe worse test MSE and
poor sparse basis recovery performance when enhancement 1 is removed. Thus, we
conclude that we need to re-privatize, _i.e.,_ add lower-order noise to the
correlations/gradients computed on the chosen feature set to obtain
competitive sparse basis recovery in the DP-FL case.

Samples | Test MSE | Number of samples recovered
---|---|---
SPriFed-OMP (low noise) |  SPriFed-OMP (high noise) |  SPriFed-OMP (low noise) |  SPriFed-OMP (high noise)
$n=2000$ | 0.83 | 335.84 | 1 | 0
$n=4000$ | 0.52 | 6.45 | 2.67 | 1.67
$n=8000$ | 0.35 | 0.77 | 7.67 | 5.67

Table 1: Enhancement 1 Performance Improvement: Comparison of Test MSE and
basis recovery performance for when Algorithm 3 runs with and without
enhancement 1. The experiments are run over 3 randomized trials over the
synthetic setup (above) with $p=10000$, $s=10$, and a privacy guarantee of
$(5.74,10^{-4})$ DP.

Importance of Enhancement 2: Privatizing the gradient

In this experiment, we look closely at the two kinds of artifacts released by
Algorithms 3 and 4, _i.e._ the correlations and gradients respectively. We
measure the change in the artifacts’ value over iterations by measuring and
plotting the total $l_{2}$ norm values for each artifact. To do so, we run a
non-private version of Algorithm 3 and plot the changing correlation norms and
gradient $l_{2}$ norms over multiple iterations. Here, in each iteration we
use the correlation value with respect to the current newly chosen feature $j$
(referred in Figure 1’s legend). We observe from Figure 1 that the gradient
norm values decrease over time while most of the correlation norms are
constant. As expected from the way OMP is designed, only the value of the
correlation norm $\lVert\bm{X_{j}^{T}y}\rVert_{2}$ decreases over iterations
since the newly chosen features have lower correlation values than previous
features. Thus, we conclude that we can employ more aggressive clipping in
Algorithm 4 than in Algorithm 3. That is, for gradients, a lower clipping
bound should work well, especially toward the final rounds, while such lower
clipping bounds will significantly affect the correlations irrespective of the
iteration.

Figure 1: Enhancement 2 Performance Intuition: Change in size of artifacts
(correlations and gradient) measured by their norm over multiple iterations.

(a) fix $s=5,n=2000$, vary $p$

(b) fix $s=5,p=2500$, vary $n$

Figure 2: Test MSE is shown for both SPriFed-OMP and DP-SGD for privacy
parameters $(4.94,10^{-4})$. Figure (a) fixes the sample size and varies the
model dimensions; Figure (b) fixes the model dimensions and varies the sample
size. Measurements averaged over $3$ randomized simulation runs.

(a) fix $s=5,p=2500$, vary $n$

(b) fix $s=5,p=10000$, vary $n$

Figure 3: The number of basis correctly recovered by SPriFed-OMP for privacy
parameters $(4.94,10^{-4})$. We choose $s=5$. Figure (a) demonstrates basis
recovery for model sparsity $p=2500$ over varying sample sizes. Figure (b)
demonstrates basis recovery for model sparsity $p=10000$ over varying sample
sizes. Measurements averaged over $3$ randomized simulation runs.

(a) fix $s=10,p=2500$, vary $n$

(b) fix $s=10,p=10000$, vary $n$

(c) fix $s=10$, vary $p$ and $n$

Figure 4: The number of basis correctly recovered by SPriFed-OMP for privacy
parameters $(5.34,10^{-4})$. We choose $s=10$, and figures (a) and (b) have an
additive error with standard deviation $\sigma_{\varepsilon}=0.001$. Figure
(a) demonstrates basis recovery for model sparsity $p=2500$ over varying
sample sizes. Figure (b) demonstrates basis recovery for model sparsity
$p=10000$ over varying sample sizes. Figure (c) varies $p=20000$ and $p=40000$
while increasing the additive error’s standard deviation to
$\sigma_{\varepsilon}=0.1$. Measurements averaged over $3$ randomized
simulation runs.

###  7.2 Performance Comparison: Synthetic Data Sets

In the rest of the section, we will compare SPriFed-OMP with various
baselines. We first discuss the experimental setup and the baselines. In
SPriFed-OMP (Algorithm 3), we set the parameter $\mu_{p}$ to be either $0.4$
(for $s=10$) or $0.543$ (for $s=5$) and set the parameter $\mu_{s}$ as $0.02$.
The cumulative privacy cost $(\epsilon,\delta)$ for all the experiments is set
to at most $(5.34,10^{-4})$ (for DP) or equivalently $1.32$ (GDP). Note that
SPriFed-OMP does not require hyperparameter tuning besides choosing the
privacy parameters $\mu_{s}$ and $\mu_{p}$. Since $p\gg s$, we can easily
tolerate $\mu_{s}$ much smaller than $\mu_{p}$. Thus, the $\mu_{p}$ parameter
will primarily dominate the privacy cost.

The first baseline DP-SGD is implemented to minimize the $L_{1}$-regularized
mean square error. After the algorithm is terminated, we pick the top-$s$
elements and leverage those to compute the final model and the test MSE. One
requirement for DP-SGD is that we also need to choose the hyper-parameters,
which are the coefficient for the $L_{1}$-regularization term and the learning
rate. Although hyper-parameter tuning may incur an additional privacy cost for
DP-SGD, here, to benefit DP-SGD, we perform the hyper-parameter tuning non-
privately using the Optuna library(Akiba et al., 2019). To maintain a
standardized notion of privacy definition, we also study the composed privacy
cost of DP-SGD using GDP. In particular, for achieving GDP in DP-SGD, we set
the privacy parameter $\mu_{DP-SGD}$ as either $0.4$ (when $s=10$) or $0.543$
(when $s=5$). $\mu_{DP-SGD}$ matches the $\mu_{p}$ parameter from SPriFed-OMP,
which is also being used to privatize $p$-dimensional vectors. DP-SGD stops
training when the overall privacy budget matches that of SPriFed-
OMP.444Although $\mu$-GDP alone provides an intuitive privacy definition, for
completeness, we also convert the final $\mu$-GDP privacy parameter to the
approximate DP parameters using the AutoDP Python library (Wang, 2023).

The second baseline DP-GCD is based on a modified DP-FL compatible version of
the algorithm presented in Mangold et al. (2023). Algorithm 1 in Mangold et
al. (2023) primarily follows from the Gradient Coordinate Descent (GCD)
algorithm. Like DP-SGD, although DP-GCD does not provide performance
guarantees on sparse basis recovery, it remains a natural algorithm when used
along with $l_{1}$-regularization objectives. However, DP-GCD in Mangold et
al. (2023) considered DP only in the server setting but not the DP-FL setting
as in this paper. Thus, we make an effort to compare with it by first
converting it to a DP-FL compatible version. We first present our converted
DP-GCD algorithm (Algorithm 5). Here, in each iteration, the maximum gradient
is picked privately in lines 5-8, and the model takes a step in only the
direction of this chosen gradient element in lines 10-11. Such gradient-based
optimization steps are repeated until the model converges or the privacy
budget is exhausted.

In our modified version presented in Algorithm 5, we also include two
enhancements to the DP-GCD algorithm from Mangold et al. (2023). First, we
convert DP-GCD to be DP-FL compatible by leveraging the NOISY-SMPC mechanism
while identifying the feature with the maximum gradient (lines 5-8). We
further ensure that much lower noise is added while following the gradient-
based approach in DP-GCD by separately re-privatizing (in line 10) the feature
gradient chosen via the maximum gradient mechanism. Thus, in summary, our
proposed modification of DP-GCD (Algorithm 5) significantly enhances the
private version in Mangold et al. (2023) by making it DP-FL compatible while
also ensuring that it has traits suitable for sparse basis recovery. Similar
to DP-SGD, we again perform hyper-parameter tuning non-privately using the
Optuna library (Akiba et al., 2019). For both DP-SGD and DP-GCD, we do not
include the cost of tuning these hyperparameters in the privacy cost, and
hence, the results we report below are more optimistic for these two
algorithms. Below we discuss our empirical findings.

Algorithm 5 Enhanced DP-GCD (DP-FL)

1:Input Parameters: Data:
$\bm{X=\\{X_{1},X_{2},...,X_{p}\\}}\in\mathbb{R}^{n\times p}$ where $X_{i}$ is
the $i^{th}$ feature, Response: $\bm{y}\in\mathbb{R}^{n}$, Learning Step Size:
$\gamma$, Total Steps: $T$, Privacy Parameters:
$\sigma_{1}=\frac{1}{\mu_{p}}$, Initialize Model: $\tilde{\alpha}=\bar{0}$,
Initialize Basis: $S=\emptyset$, Gradient Bound: $C$ (true bound over all
gradients or bound achieved via gradient clipping).

2:Output: Predicted Sparse Model $\tilde{\alpha}$, Predicted Sparse Basis $\S$

3:procedure DP-GCD($\bm{X},\bm{y},n,T$)

4: for $t=1:T$ do

5: for $j=0:p$ do

6: $g_{j}\leftarrow\text{{NOISY-
SMPC}}(X_{j},y-X_{S_{\mathcal{A}}}\tilde{\alpha},n,\sigma_{1}C\sqrt{p})$

7: end for

8: $j^{*}\leftarrow\operatorname*{arg\,max}_{j}|g_{j}|$

9: $S\leftarrow S\cup\\{j^{*}\\}$

10: $g_{j^{*}}^{{}^{\prime}}\leftarrow\text{{NOISY-
SMPC}}(X_{j^{*}},y-X_{S_{\mathcal{A}}}\tilde{\alpha},n,\sigma_{1}C)$

11: $\tilde{\alpha}[j^{*}]=\tilde{\alpha}[j^{*}]-\gamma
g_{j^{*}}^{{}^{\prime}}$

12: end for

13: Return $\tilde{\alpha},S$

14:end procedure

Fix $n$, Change $p$: In Fig. 2(a), we fix $s=5$ and $n=2000$ and vary
$p\in\\{2500,4000,6000,8000,10000\\}$. With rising $p$, we observe that the
test MSE loss increases much more dramatically for DP-SGD (note the
logarithmic y-axis) than SPriFed-OMP. In fact, SPriFed-OMP’s performance is
almost constant as $p$ varies since its empirical risk (Theorem 10) is mostly
dependent on $s$. Further, although not shown in the figure, SPriFed-OMP
extracts most of the basis elements $(\geq 3)$ even when $p=10000$. Finally,
even though the sample size can be quite small in comparison to the model
dimension, we obtain highly accurate models for SPriFed-OMP with test MSE
close to zero, which is in line with our theory (see Theorems 9 and 10).
Similarly, even though DP-GCD performs better than DP-SGD in terms of test MSE
(thanks to the lower noise that our enhanced DP-GCD algorithm adds), its test
MSE is still much worse than SPriFed-OMP. We will discuss this limitation of
DP-SGD in the latter part of this section.

Fix $p$, Change $n$: In Fig. 2(b), we set $s=5$ and $p=2500$ and vary
$n\in\\{400,800,1200,1600,2000\\}$. We observe that when $p\gg n$, the
performance of all algorithms is impacted. However, increasing the sample size
improves performance more dramatically for SPriFed-OMP and SPriFed-OMP-GRAD
(again, note the logarithmic scale on the y-axis). Notably, SPriFed-OMP and
SPriFed-OMP-GRAD significantly outperform DP-SGD and DP-GCD for all choices of
$n$. We specifically note that the test loss for both versions of SPriFed-OMP
gets very close to zero (as $n$ increases), which matches our theoretical
results (Theorems 9 and 10).

Fix $s$, Change $p$ and $n$: In Figs. 3 and 4, we report the number of non-
zero basis in $\alpha_{*}$ that are correctly recovered by the proposed
algorithms. We show two plots, Fig. 3 for $s=5$ and Fig. 4 for $s=10$. The
sample sizes are varied over $\\{500,1000,2000,4000,8000\\}$ and the model
dimensions are varied over $\\{2500,100000\\}$ for both $s=5$ and $s=10$.
Additionally, we also vary p over $\\{20000,400000\\}$ for the $s=10$ case.
First, we can see for correctly recovering majority of the sparse basis, the
value of $n$ increases much slower than $p$. Let us examine more closely the
cases when we vary $s$ as either $5$ or $10$. We first look at the case when
$s=5$. For $p=2500$ by setting $n=2000$, we manage to recover $3+$ out of $5$
basis elements with SPriFed-OMP-Grad. It turns out that even when $p$ is
increased to 10000, using $n=2000$ is still sufficient to recover 3 out of 5
basis. When $s=10$, we consider the number of samples required to recover $7+$
out of $10$ basis elements. For $p=2500$, $n=2000$ suffices. But when $p$ is
increased $4$ times to $10000$, we only need $n$ to be doubled to $4000$.
Finally, when $p$ is either $20000$ or $40000$, we only need $n=8000$ samples
to recover the same number of basis elements. In summary, we find that, as the
model dimension and the model sparsity increase, the sample size required for
sparse recovery increases as well, which is consistent with our main results,
Theorems 7 and 8. More interestingly, the value of $n$ increases much slower
than that of $p$. This aligns with our theoretical result that $n$ only needs
to be order-$\sqrt{p}$.

Second, as seen while comparing Figures 4(a-b) with Figure 4(c), increasing
the additive error variability affects SPriFed-OMP the most, while both the
gradient-based methods SPriFed-OMP-Grad and DP-GCD are minimally affected.
Note that the results shown in Figures 4(a-b) have system error standard
deviation of $0.001$. Here, we can see that the performance of SPriFed-OMP is
quite comparable to both SPriFed-OMP-GRAD and DP-GCD. SPriFed-OMP even manages
to outperform DP-GCD for the higher sample size range. However, Figure 4(c)
uses a standard deviation of $0.1$. Here, the performance of SPriFed-OMP dips
significantly in comparison to both SPriFed-OMP-GRAD and DP-GCD. Thus, error
variability affects SPriFed-OMP more than the other methods. These results are
in line with our theory (refer Theorems 7 and 8) that shows that SPriFed-OMP
has stricter system error requirements compared to SPriFed-OMP-GRAD.

Let us now discuss the sparse basis recovery performance of the two baselines.
Note that both DP-GCD and DP-SGD converge to a private model by leveraging
gradient-based optimization. However, DP-GCD adds much lower noise than DP-SGD
since, at each round, it only updates a single dimension rather than all
dimensions as done by DP-SGD. This explains why we see a better performance
for DP-GCD than DP-SGD. However, DP-GCD still adds significant noise to the
model compared to both SPriFed-OMP and SPriFed-OMP-GRAD due to the following
two reasons. First, due to its gradient-based optimization, DP-GCD may take
more than one step to optimize any given direction fully. In contrast, for
each new basis, SPriFed-OMP only collects new correlation values (or
gradients) once. Thus, depending on the number of steps taken, the amount of
noise added by DP-GCD per dimension is often much higher. For instance, if the
same dimension is picked $t$ times, then the noise added will be of $t$ orders
higher for DP-GCD than the SPriFed-OMP class of methods. The remark in this
section covers further details. Second, unlike SPriFed-OMP and SPriFed-OMP-
GRAD, DP-GCD does not compute the exact true value of the model (using the OLS
formula) in each round rather than choosing to add noise to the previously
predicted model continually and thus the noise added in previous rounds might
affect the predicted model significantly. In contrast, private OMP estimates
the model parameters $\alpha$ directly using the OLS formula (Algorithm 2)
based on the re-privatized correlation values in Lines (14)) and (9) of
Algorithms 3 and 4 respectively, which have much lower noise. Thus, this noise
accumulation effect is greatly reduced.

Remark (Significant difference of Test MSE between DP-SGD, DP-GCD and both
flavors of SPriFed-OMP): From Figures 2(a) and 2(b), we observe that the test
MSE of SPriFed-OMP and DP-SGD are of different orders (note the logarithmic
scale on the y-axis). This difference is because, in SPriFed-OMP, we only
introduce noise with a variance of order-$s$ to the final model. In contrast,
DP-SGD adds noise of order-$p$ variance to all its model dimensions.
Therefore, even though DP-SGD selects a subset of these dimensions for the
final model, the considerable noise significantly impacts the resulting model
coefficients. As a result, the test MSE for DP-SGD is much worse than the one
we obtain from SPriFed-OMP. Once more, we note that due to such a high noise
variance, DP-SGD cannot recover any of the true basis in our experiments. For
the case of DP-GCD, the noise added to the model is of order $s$ as well.
However, DP-GCD follows a gradient-based iterative approach where noise is
added to the gradient. Then, this noise is passed to the model based on the
feature selected in that step, the learning rate, and the number of total
steps taken. Unlike DP-GCD, the SPriFed-OMP flavors add noise to the exact
value of the true model based on the currently selected features. Thus,
overall, SPriFed-OMP only adds noise via the privacy mechanism. On the other
hand, DP-GCD adds a similar noise via the privacy mechanism while also
including a bias term based on the difference between the DP-GCD model and the
optimal linear model for the chosen feature set.

###  7.3 Performance Comparison: More Realistic Data Sets

In the second set of experiments, we test our proposed algorithms on the
synthetic MADELON dataset from the NIPS 2003 challenge (Guyon et al., 2004), a
popular dataset used in the sparse recovery literature. The original MADELON
dataset has $20$ true features, but only $5$ of these are true raw features,
and the remaining $15$ are linear combinations of these raw features. Its
total model dimension is $p=500$, with $n=2000$ training samples. However,
since we are only interested in the $p>n$ case, we randomly sub-sample $300$
samples in each trial. We run $5$ randomized trials and report the average
results. Here, we set $\mu_{p}=0.45,\mu_{s}=0.09$ for SPriFed-OMP and
$\mu_{DP-SGD}=0.45$ for DP-SGD. The corresponding DP guarantee is
$(8.38,10^{-3})$, and the GDP guarantee is $1.72$. For SPriFed-OMP, we clip
the elements of $\bm{X}$ to $0.12$ and elements of $\bm{y}$ to $0.36$. We clip
the gradient elements to $1$. As expected, DP-SGD performs the worst with a
high Test MSE of 105.73. DP-GCD performs reasonably well with Test MSE of
1.23. Finally, SPriFed-OMP has a Test MSE of 1.55 while SPriFed-OMP-GRAD
handily outperforms DP-GCD with a Test MSE of only 0.27. We attribute the
relatively better performance of DP-GCD for this data set to two reasons, (1)
proper hyperparameter tuning of DP-GCD; thus, we expect that the performance
of DP-GCD would be worse if the privacy cost of hyperparameter tuning must
also be accounted for; and (2) due to the synthetic nature of MADELON, DP-GCD
might perform well even when it does not pick the exact best features (since
MADELON has 15 features that are linear combinations of the true 5 features).

In the third set of experiments, we show results on two real high-dimensional
datasets from the repository presented in Drysdale (2022). Here, we consider
the high-dimensional datasets chop and gse1992. The chop dataset contains
$414$ total samples with $3836$ features, and gse1992 contains $124$ samples
with $15537$ features. Similar to MADELON, we sub-sample the dimensions of
chop to be $2000$ and that of gse1992 to be $500$ in each trial. The rest of
the privacy setup is identical to MADELON. We present the results for these
datasets in Table 2. We note that both SPriFed-OMP and SPriFed-OMP-GRAD
outperforms all other algorithms over both datasets with SPriFed-OMP-GRAD
providing the best utility. As previously noted in the remark above, we
attribute the success of both versions of SPriFed-OMP to their ability to add
significantly lower order noise to the true model value while leveraging the
true value of the model using the ordinary least squares formula. Note that
DP-GCD also adds lower-order noise to the model. However, unlike both flavors
of SPriFed-OMP, DP-GCD does not add noise to the exact value of the true
model, instead relying on a gradient-based approach that iteratively adds
noise to the gradient. The gradient-based noise is then passed to the model
value based on the feature selected and the learning rate.

Algorithm | chop | gse1992
---|---|---
SPriFed-OMP | 0.101 | 1.90
SPriFed-OMP-GRAD | 0.154 | 1.15
DP-GCD | 18.79 | 525.60
DP-SGD | 546.76 | 988.72
Table 2: Test MSE Performance: We observe that without model/basis selection,
both DP-SGD and DP-GCD are required to add substantial noise in the models,
leading to much higher Test MSE than both SPriFed-OMP and SPriFed-OMP-GRAD.
However, both SPriFed-OMP and SPriFed-OMP-GRAD enjoy close to optimal test
MSE.

##  8 Conclusion

In this work, we propose two new private sparse basis recovery algorithms
called SPriFed-OMP and SPriFed-OMP-GRAD for the DP-FL setting based on the
Orthogonal Matching Pursuit (OMP) algorithm. Specifically, we prove
analytically that both algorithms can recover the sparse basis in a finite
number of steps. Further, we bound the privacy cost while leveraging the
obtained low-dimensional model to obtain dimension-free model estimation error
and empirical risk. Here, SPriFed-OMP is our first attempt at leveraging OMP
for sparse basis recovery in the DP-FL setting. In contrast, SPriFed-OMP-GRAD
improves upon SPriFed-OMP (both analytically and empirically) by requiring a
simpler DP-FL mechanism, fewer samples, and a less complicated analysis. To
the best of our knowledge, these are the first DP-FL algorithms that
successfully manage sparse recovery under the RIP assumption while only
requiring $n=\mathcal{O}(\sqrt{p})$ samples. Our experimental results on both
synthetic and real datasets strongly corroborate our theoretical analysis. For
future work, we will explore DP-FL algorithms that can perform sparse recovery
even without the RIP assumption.

## References

  * Abadi et al. (2016) Martin Abadi, Andy Chu, Ian Goodfellow, H Brendan McMahan, Ilya Mironov, Kunal Talwar, and Li Zhang.  Deep learning with differential privacy.  In _Proceedings of the 2016 ACM SIGSAC conference on computer and communications security_ , pp. 308–318, 2016.
  * Akiba et al. (2019) Takuya Akiba, Shotaro Sano, Toshihiko Yanase, Takeru Ohta, and Masanori Koyama.  Optuna: A next-generation hyperparameter optimization framework.  In _Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining_, pp. 2623–2631, 2019.
  * Asi et al. (2021) Hilal Asi, Vitaly Feldman, Tomer Koren, and Kunal Talwar.  Private stochastic convex optimization: Optimal rates in l1 geometry.  In _International Conference on Machine Learning_ , pp. 393–403. PMLR, 2021.
  * Bassily et al. (2014) Raef Bassily, Adam Smith, and Abhradeep Thakurta.  Private empirical risk minimization: Efficient algorithms and tight error bounds.  In _2014 IEEE 55th annual symposium on foundations of computer science_ , pp. 464–473. IEEE, 2014.
  * Bassily et al. (2021) Raef Bassily, Cristóbal Guzmán, and Anupama Nandi.  Non-euclidean differentially private stochastic convex optimization.  In _Conference on Learning Theory_ , pp. 474–499. PMLR, 2021.
  * Belloni et al. (2014) Alexandre Belloni, Victor Chernozhukov, and Christian Hansen.  High-dimensional methods and inference on structural and treatment effects.  _Journal of Economic Perspectives_ , 28(2):29–50, 2014.
  * Bishop et al. (1995) Christopher M Bishop et al.  _Neural networks for pattern recognition_.  Oxford university press, 1995.
  * Bonawitz et al. (2017) Keith Bonawitz, Vladimir Ivanov, Ben Kreuter, Antonio Marcedone, H Brendan McMahan, Sarvar Patel, Daniel Ramage, Aaron Segal, and Karn Seth.  Practical secure aggregation for privacy-preserving machine learning.  In _proceedings of the 2017 ACM SIGSAC Conference on Computer and Communications Security_ , pp. 1175–1191, 2017.
  * Boucheron et al. (2013) Stéphane Boucheron, Gábor Lugosi, and Pascal Massart.  _Concentration inequalities: A nonasymptotic theory of independence_.  Oxford university press, 2013.
  * Bun & Steinke (2016) Mark Bun and Thomas Steinke.  Concentrated differential privacy: Simplifications, extensions, and lower bounds.  In _Theory of Cryptography Conference_ , pp. 635–658. Springer, 2016.
  * Candes & Tao (2005) Emmanuel J Candes and Terence Tao.  Decoding by linear programming.  _IEEE transactions on information theory_ , 51(12):4203–4215, 2005.
  * Chen & Caramanis (2012) Yudong Chen and Constantine Caramanis.  Orthogonal matching pursuit with noisy and missing data: Low and high dimensional results.  _arXiv preprint arXiv:1206.0823_ , 2012.
  * Clarke et al. (2008) Robert Clarke, Habtom W Ressom, Antai Wang, Jianhua Xuan, Minetta C Liu, Edmund A Gehan, and Yue Wang.  The properties of high-dimensional data spaces: implications for exploring gene and protein expression data.  _Nature reviews cancer_ , 8(1):37–49, 2008.
  * Dong et al. (2019) Jinshuo Dong, Aaron Roth, and Weijie J Su.  Gaussian differential privacy.  _arXiv preprint arXiv:1905.02383_ , 2019.
  * Drysdale (2022) Erik Drysdale.  SurvSet: An open-source time-to-event dataset repository.  _arXiv preprint arXiv:2203.03094_ , 2022.
  * Dwork et al. (2014) Cynthia Dwork, Aaron Roth, et al.  The algorithmic foundations of differential privacy.  _Found. Trends Theor. Comput. Sci._ , 9(3-4):211–407, 2014.
  * Efron et al. (2004) Bradley Efron, Trevor Hastie, Iain Johnstone, and Robert Tibshirani.  Least angle regression.  2004\.
  * Fan et al. (2011) Jianqing Fan, Jinchi Lv, and Lei Qi.  Sparse high-dimensional models in economics.  _Annu. Rev. Econ._ , 3(1):291–317, 2011.
  * Geiping et al. (2020) Jonas Geiping, Hartmut Bauermeister, Hannah Dröge, and Michael Moeller.  Inverting gradients-how easy is it to break privacy in federated learning?  _Advances in Neural Information Processing Systems_ , 33:16937–16947, 2020.
  * Guyon et al. (2004) Isabelle Guyon, Steve Gunn, Asa Ben-Hur, and Gideon Dror.  Result analysis of the nips 2003 feature selection challenge.  _Advances in neural information processing systems_ , 17, 2004.
  * Hu et al. (2022) Lijie Hu, Shuo Ni, Hanshen Xiao, and Di Wang.  High dimensional differentially private stochastic optimization with heavy-tailed data.  In _Proceedings of the 41st ACM SIGMOD-SIGACT-SIGAI Symposium on Principles of Database Systems_ , pp. 227–236, 2022.
  * Huang et al. (2021) Yangsibo Huang, Samyak Gupta, Zhao Song, Kai Li, and Sanjeev Arora.  Evaluating gradient inversion attacks and defenses in federated learning.  _Advances in Neural Information Processing Systems_ , 34:7232–7241, 2021.
  * Jalali et al. (2011) Ali Jalali, Christopher Johnson, and Pradeep Ravikumar.  On learning discrete graphical models using greedy methods.  _Advances in Neural Information Processing Systems_ , 24, 2011.
  * Kairouz et al. (2015) Peter Kairouz, Sewoong Oh, and Pramod Viswanath.  The composition theorem for differential privacy.  In _International conference on machine learning_ , pp. 1376–1385. PMLR, 2015.
  * Kairouz et al. (2021) Peter Kairouz, Ziyu Liu, and Thomas Steinke.  The distributed discrete gaussian mechanism for federated learning with secure aggregation.  In _International Conference on Machine Learning_ , pp. 5201–5212. PMLR, 2021.
  * Kifer et al. (2012) Daniel Kifer, Adam Smith, and Abhradeep Thakurta.  Private convex empirical risk minimization and high-dimensional regression.  In _Conference on Learning Theory_ , pp. 25–1. JMLR Workshop and Conference Proceedings, 2012.
  * Liang & Kelemen (2008) Yulan Liang and Arpad Kelemen.  Statistical advances and challenges for analyzing correlated high dimensional snp data in genomic study for complex diseases.  2008\.
  * Mangold et al. (2023) Paul Mangold, Aurélien Bellet, Joseph Salmon, and Marc Tommasi.  High-dimensional private empirical risk minimization by greedy coordinate descent.  In _International Conference on Artificial Intelligence and Statistics_ , pp. 4894–4916. PMLR, 2023.
  * McMahan et al. (2017) Brendan McMahan, Eider Moore, Daniel Ramage, Seth Hampson, and Blaise Aguera y Arcas.  Communication-efficient learning of deep networks from decentralized data.  In _Artificial intelligence and statistics_ , pp. 1273–1282. PMLR, 2017.
  * Meinshausen (2007) Nicolai Meinshausen.  Relaxed lasso.  _Computational Statistics & Data Analysis_, 52(1):374–393, 2007.
  * Meinshausen & Yu (2009) Nicolai Meinshausen and Bin Yu.  Lasso-type recovery of sparse representations for high-dimensional data.  2009\.
  * Mironov (2017) Ilya Mironov.  Rényi differential privacy.  In _2017 IEEE 30th computer security foundations symposium (CSF)_ , pp. 263–275. IEEE, 2017.
  * Needell & 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.
  * Petersen et al. (2008) Kaare Brandt Petersen, Michael Syskind Pedersen, et al.  The matrix cookbook.  _Technical University of Denmark_ , 7(15):510, 2008.
  * Rudelson & Vershynin (2008) Mark Rudelson and Roman Vershynin.  The littlewood–offord problem and invertibility of random matrices.  _Advances in Mathematics_ , 218(2):600–633, 2008.
  * Talwar et al. (2015) Kunal Talwar, Abhradeep Guha Thakurta, and Li Zhang.  Nearly optimal private lasso.  _Advances in Neural Information Processing Systems_ , 28, 2015.
  * Thakurta & Smith (2013) Abhradeep Guha Thakurta and Adam Smith.  Differentially private feature selection via stability arguments, and the robustness of the lasso.  In _Conference on Learning Theory_ , pp. 819–850. PMLR, 2013.
  * Tropp et al. (2007) Joel Tropp, Anna C Gilbert, et al.  Signal recovery from partial information via orthogonal matching pursuit.  _IEEE Trans. Inform. Theory_ , 53(12):4655–4666, 2007.
  * Vershynin (2010) Roman Vershynin.  Introduction to the non-asymptotic analysis of random matrices.  _arXiv preprint arXiv:1011.3027_ , 2010.
  * Vershynin (2018) Roman Vershynin.  _High-dimensional probability: An introduction with applications in data science_ , volume 47.  Cambridge university press, 2018.
  * Wang & Xu (2019) Di Wang and Jinhui Xu.  On sparse linear regression in the local differential privacy model.  In _International Conference on Machine Learning_ , pp. 6628–6637. PMLR, 2019.
  * Wang & Gu (2019) Lingxiao Wang and Quanquan Gu.  Differentially private iterative gradient hard thresholding for sparse learning.  In _28th International Joint Conference on Artificial Intelligence_ , 2019.
  * Wang (2023) Yu-Xiang Wang.  autodp: A flexible and easy-to-use package for differential privacy.  https://github.com/yuxiangw/autodp, 2023.
  * Wasserman & Roeder (2009) Larry Wasserman and Kathryn Roeder.  High dimensional variable selection.  _Annals of statistics_ , 37(5A):2178, 2009.
  * Zhang (2009) Tong Zhang.  On the consistency of feature selection using greedy least squares regression.  _Journal of Machine Learning Research_ , 10(3), 2009.
  * Zhang (2011) Tong Zhang.  Adaptive forward-backward greedy algorithm for learning sparse representations.  _IEEE transactions on information theory_ , 57(7):4689–4708, 2011.
  * Zhao & Yu (2006) Peng Zhao and Bin Yu.  On model selection consistency of lasso.  _The Journal of Machine Learning Research_ , 7:2541–2563, 2006.

##  Appendix A Differential Privacy

Differential privacy was introduced to privatize changes caused by individual
clients (Dwork et al., 2014). For a given function, differential privacy
obscures the maximum possible change in the function’s output due to the
inclusion or removal of a single client’s data. Here, obscuring is allowed by
introducing noise and, thus, a form of deniability in the function’s output.
Deniability allows an individual client to deny its presence in the dataset
and disallows an adversary from linking the client or its data to a non-
privatized record. Differential privacy is significantly stronger than
previous anonymization methods as it works even in the worst-case scenario
when the entire dataset is publicly released except for a single client’s
data. Here, we look at the precise mathematical definition of differential
privacy, simple techniques to introduce privacy, and its additional
properties.

######  Definition A.1.

[Approximate Differential Privacy] For any two neighboring datasets
$X,X{{}^{\prime}}\subseteq\mathcal{S}$ differing over a single sample, we say
that the randomized mechanism $\mathcal{M}:\mathcal{S}\rightarrow\mathcal{R}$
is $(\epsilon,\delta)$ differentially private if,

$\displaystyle\text{Pr}[\mathcal{M}(X)\in\mathcal{R}]\leq
e^{\epsilon}\cdot\text{Pr}[\mathcal{M}(X^{{}^{\prime}})\in\mathcal{R}]+\delta$

######  Definition A.2 ($L_{2}$ sensitivity of a function).

Given a deterministic function $f$, consider all possible pairs of conformable
neighboring datasets in the domain of $f$ differing on a single row. Then we
denote $f$’s $L_{2}$ sensitivity as

$\displaystyle\Delta_{2}f\coloneqq\max_{x,x\in
dom(f),||x}||f(x)-f(x^{\prime})||_{2}$

Thus, $\Delta_{2}f$ records the maximum possible change in $f$ due to a change
in a single row. Often each row belongs to a different client and thus,
$\Delta_{2}f$ will record the maximum change due to a change in a client or
the inclusion or removal of a client from the dataset.

We can achieve differential privacy by the Laplacian, Exponential, and
Gaussian mechanisms (Dwork et al., 2014). Here, we use the Gaussian mechanism
as it performs better for multi-step (or composed) privacy algorithms. In
particular, we leverage the Gaussian Differential Privacy (GDP) composition
method for optimally tracking the privacy loss over multiple steps (Dong et
al., 2019). Although, other composition methods such as advanced composition
(Kairouz et al., 2015), Renyi DP (Mironov, 2017), and concentrated DP (Bun &
Steinke, 2016) exist, we chose GDP for its superior composition performance
for Gaussian mechanisms. For completeness, we state the connection between the
Gaussian mechanism and the GDP privacy constant.

######  Lemma 11.

[GDP Gaussian mechanism (Theorem 2.7 from (Dong et al., 2019))] Let statistic
$f(\cdot)$ be computed over the dataset $\bm{X}$. Then, the randomized
Gaussian mechanism $\mathcal{M}(\bm{X})=f(\bm{X})+\eta$ where
$\eta\sim\mathcal{N}(0,\frac{\Delta_{2}f^{2}}{\mu^{2}})$ is $\mu-GDP$.

######  Lemma 12.

[Restated Theorem 2.7 from Dong et al. (2019)] Let statistic $f(\cdot)$ be
computed over the dataset $\bm{X}$. Then, the randomized Gaussian mechanism
$\mathcal{M}(\bm{X})=f(\bm{X})+\eta$ where
$\eta\sim\mathcal{N}(0,\frac{\Delta_{2}f^{2}}{\mu^{2}})$ is $\mu-GDP$.

An important property of DP, known as post-processing maintains that no amount
of processing on a differentially private system will cause the new output to
leak any more information than the unprocessed input. This is expressed
clearly in Lemma 13.

######  Lemma 13 (Differential Privacy Post-Processing).

Consider a randomized mechanism $M$ that is $(\varepsilon,\delta)$
differentially private. Now, let $M_{post}$ be another randomized mechanism
that can be applied to the output of $M$. Then, we state that the output of
the combination of $M_{post}(M(\cdot))$ function group has at least
differential privacy guarantee of $(\varepsilon,\delta)$ or higher.

##  Appendix B Compressed Sensing Theory

Let $\Omega=\\{1,2,...,p\\}$ be the overall feature index set and the support
$S\subseteq\Omega$ of $\beta^{*}$ be defined such that,
$S=\\{i|i\in\Omega,\beta^{*}_{i}\neq 0\\}$. By definition, $|S|\leq s$ (where
$s$ is the sparsity index of the dataset). For a vector $v\in\mathbb{R}^{p}$,
let us define $v_{W}$ as the restriction of $v$ to the indices in $W$ where
$W\subseteq\Omega$. Similarly, for a matrix $Q\in\mathbb{R}^{n\times p}$ let
us define the sub-matrix $Q_{W}\in\mathbb{R}^{n\times|W|}$ s.t. it only
contains columns indexed by $W(\subseteq\Omega)$. Our primary assumption for
the data matrix follows the Restricted Isometry Property (RIP) as detailed in
Assumption 1.

Intuitively, assumption 1 for $X$, a matrix satisfying RIP of order $K$ states
that the eigenvalues for all $K$-dimensional (or smaller) sub-matrices derived
from $X$, lie in a tight interval around $1$. Interestingly, the assumption
holds for several popular sub-Gaussian random matrices. Construction of such
random RIP matrices is detailed in Theorem 5.65 from Vershynin (2010). We note
that although we use the random matrix construction as a reference for our
primary analysis, we only require a matrix (even deterministic) normalized by
$\sqrt{n}$ that satisfies the RIP condition and and it should contain
(probably) bounded values. Further details are omitted due to page limits.

######  Lemma 14.

[Construction of Random Matrices satisfying RIP] (Theorem 5.65 from Vershynin
(2010)) Let $A\in\mathbb{R}^{n\times p}$ be a subgaussian random matrix with
independent and isotropic rows. Then the normalized matrix
$\frac{A}{\sqrt{n}}$ satisfies the following for every sparsity level $1\leq
k\leq p$ and every number $\delta\in(0,1)$,

$\displaystyle
ifn\geq\frac{Ck\log(\frac{ep}{k})}{\delta^{2}}\rightarrow\delta_{k}(\frac{A}{\sqrt{n}})\leq\delta$

with a probability at least $1-2e^{-c\delta^{2}n}$. Here, $C=C_{K},c=c_{K}>0$
depend only on the subgaussian norm $K=\max_{i}||A_{i}||_{\Psi_{2}}$ of the
rows of $A$. Here, $k$ is the RIP order, and $\delta_{k}(\cdot)$ is the RIC.
The exact construction is demonstrated below.

One instance of such a random matrix satisfying RIP, $X\in\mathbb{R}^{n\times
p}$, can be constructed with the following steps,

  1. 1.

For a given row $v$ in $X$, independently sample the row’s elements from a
sub-Gaussian distribution s.t. each row is isotropic. $v\in\mathbb{R}^{p}$ is
an isotropic vector if $E[vv^{T}]=\mathbb{I}$.

  2. 2.

The standard normal distribution satisfies the isotropic property. Therefore,
assume $X_{ij}\sim\mathcal{N}(0,1),\text{for all }i\in[n],j\in[p]$. We note
that alternative distributions as detailed in Vershynin (2010) also satisfy
the isotropy property.

  3. 3.

Normalize $X\rightarrow\frac{X}{\sqrt{n}}$.

Given an RIP satisfying matrix $X\in\mathbb{R}^{n\times p}$ satisfying RIP of
order $s$ with an RIC of constant $\zeta_{s}$, we note some useful properties
(Lemmas 15-18) derived from the RIP condition ((Candes & Tao, 2005)).

Consider the set of vectors $\mathcal{W}\coloneqq\\{w\in\mathbb{R}^{p}$ | $||w||_{0}\leq s\\}$. The support for vectors in $\mathcal{W}$ is $\mathcal{S}$. For a vector $w\in\mathcal{W}$ and by RIP and using $||X_{S}w_{S}||_{2}=||Xw||2$,

$\displaystyle(1-\zeta_{s})||w||_{2}^{2}$
$\displaystyle\leq||X_{S}w_{S}||_{2}^{2}\leq(1+\zeta_{s})||w||_{2}^{2}$
$\displaystyle\rightarrow(1-\zeta_{s})||w||_{2}^{2}$
$\displaystyle\leq||Xw||_{2}^{2}\leq(1+\zeta_{s})||w||_{2}^{2}$
$\displaystyle\rightarrow(1-\zeta_{s})$
$\displaystyle\leq\frac{w^{T}X^{T}Xw}{w^{T}w}\leq(1+\zeta_{s})$ (2)

Thus, the following bound immediately follows,

$\displaystyle(1-\zeta_{s})$
$\displaystyle\leq\sigma_{min}(X_{S}^{T}X_{S}),\sigma_{\max}(X_{S}^{T}X_{S})\leq(1+\zeta_{s})$

######  Lemma 15.

[Monotonocity of the Isometric Constant] If a measurement matrix
$X\in\mathbb{R}^{n\times p}$ satisfies RIP of orders $K_{1}$ and $K_{2}$ s.t.
$K_{1}\leq K_{2}$ then their corresponding RICs follow the same order _i.e.,_
$\zeta_{K_{1}}\leq\zeta_{K_{2}}$.

######  Lemma 16.

[Direct RIP Consequences] Let $S\subseteq\Omega$ be the support of the true
model $\beta^{*}$. If $\zeta_{s}<1;s=|S|$ then we can show that for the
measurement matrix $X$ satisfying RIP with order $s$ and for any
$w\in\mathbb{R}^{s}$,

$\displaystyle(1-\zeta_{s})||w||_{2}\leq||X_{S}^{T}X_{S}w||_{2}\leq(1+\zeta_{s})||w||_{2}$
$\displaystyle\frac{1}{1+\zeta_{s}}||w||_{2}\leq||(X_{S}^{T}X_{S})^{-1}w||_{2}\leq\frac{1}{1-\zeta_{s}}||w||_{2}$

######  Lemma 17.

[RIP: Near Orthogonality] Let $S_{1},S_{2}\subseteq\Omega s.t.,S_{1}\cap
S_{2}=\emptyset$. If $\zeta_{|S_{1}|+|S_{2}|}<1$ and $S_{1}\cup S_{2}=S$ then
for the measurement matrix $X$ satisfying RIP of order $s$ and for
$u\in\mathbb{R}^{|S_{2}|}$ we have

$\displaystyle||X_{S_{1}}^{T}X_{S_{2}}u||_{2}\leq\zeta_{|S_{1}|+|S_{2}|}||u||_{2}$

######  Lemma 18.

Let $\mathcal{S}\subseteq\Omega$ be the support of $\beta^{*}$ (our ground
truth model parameter). If $\zeta_{s}<1;s=|\mathcal{S}|$ then we can show for
the measurement matrix $X$ satisfying RIP with order $s$ and for any
$w\in\mathbb{R}^{n}$.

$\displaystyle||X_{\mathcal{S}}^{T}w||_{2}\leq\sqrt{1+\zeta_{s}}||w||_{2}$

##  Appendix C Useful Results for Privacy Analysis

Proof of Lemma 4 Consider two neighboring input-output pairs $X,y$ and
$X^{{}^{\prime}},y^{{}^{\prime}}$ such that the pairs differ in a single
sample at row $k$. We approach this proof by first considering the sensitivity
of the $j^{th}$ column ($X_{j}$) and then computing the combined sensitivity
of the $p$ columns.

$\displaystyle\Delta_{2}(X_{j}y_{j})$
$\displaystyle=\max_{k\in[n]}||(X_{j}^{{}^{\prime}})^{T}y^{{}^{\prime}}-X_{j}^{T}y||_{2}$
$\displaystyle\leq 2\max_{k\in[n]}||X_{jk}y_{k}||_{2}$ $\displaystyle\leq
2\max_{k\in[n]}X_{jk}y_{j}\leq X_{M}y_{M}$
$\displaystyle\rightarrow\Delta_{2}(X^{T}y)$ $\displaystyle\leq
2\sqrt{p}X_{M}y_{M}=\mathcal{O}(\sqrt{p})$

Similarly, let us consider the sensitivity of $X^{T}X_{S_{\mathcal{A}}}$,

$\displaystyle\Delta_{2}(X^{T}X_{j})$
$\displaystyle=\max_{k\in[n]}||(X^{{}^{\prime}})^{T}X^{{}^{\prime}}_{i}-X^{T}X_{i}||_{2}$
$\displaystyle\leq 2\max_{k\in[n]}(||[X_{lk}X_{jk}]_{l\in[p]}||_{2})$
$\displaystyle\leq 2\sqrt{p}X_{M}^{2}$
$\displaystyle\rightarrow\Delta_{2}(X^{T}X_{S_{\mathcal{A}}})$
$\displaystyle\leq 2X_{M}^{2}\sqrt{ps}=\mathcal{O}(\sqrt{ps})$

##  Appendix D Useful Results for optimality of SPriFed-OMP

######  Lemma 19.

Consider the term $M=[\mathbb{I}+NA^{-1}]^{-1}v$ where
$N\coloneqq\frac{\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}}{n},A\coloneqq(\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n})^{-1}$,
$v$ is any conformable vector and $n$ is the sample size. The remaining terms
have been borrowed from Algorithm 3. Note that each element in
$\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}$ has a distribution of
$\mathcal{N}(0,p/\mu_{p}^{2})$. We show that

$\displaystyle||M||_{2}\leq\kappa_{M}||v||_{2};\kappa_{M}$
$\displaystyle=\frac{1}{1+\frac{\varepsilon_{\eta}\sqrt{s}\kappa_{s}}{n(1-\zeta_{s+1})}}$

Here, $\varepsilon_{\eta}$ is chosen as per the lower bound of the smallest
singular value of the random matrix $X_{S_{\mathcal{A}}}$ (by Theorem 1.2 from
Rudelson & Vershynin (2008)), $n=\sqrt{ps}\kappa_{n}$ from assumptions in
Theorem 7 and so we have that each row in
$\frac{\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}}{n}$ has the
distribution $\sim\frac{1}{(n/\sqrt{s}\kappa_{s})}\mathcal{N}(0,1)$.

Proof:

$\displaystyle||M||_{2}$ $\displaystyle=||[\mathbb{I}+NA^{-1}]^{-1}v||_{2}$
$\displaystyle\overset{(1)}{\leq}\frac{||v||_{2}}{\sigma_{\min}[\mathbb{I}+NA^{-1}]}$
$\displaystyle\overset{(2)}{=}\frac{||v||_{2}}{1+\sigma_{\min}[NA^{-1}]}$
$\displaystyle\overset{(3)}{\leq}\frac{||v||_{2}}{1+\sigma_{\min}[N]\sigma_{\min}[A^{-1}]}$
$\displaystyle\overset{(4)}{\leq}\frac{||v||_{2}}{1+\frac{s\varepsilon_{\eta}\kappa_{s}}{n\sqrt{s}}(1/\sigma_{\max}[A])}$
$\displaystyle\overset{(5)}{\leq}\frac{||v||_{2}}{1+\frac{\sqrt{s}\varepsilon_{\eta}\kappa_{s}}{n(1+\zeta_{s+1})}}$

Here (1) follows from Lemma 30 and since
$\sigma_{\max}[K^{-1}]=\frac{1}{\sigma_{\min}[K]}$. (2) is true since
$||k\mathbb{I}+K||_{2}=k+||K||_{2}$ where $\mathbb{I}$ is the identity matrix
and $K$ is some conformable matrix. Furthermore, with high probability, the
minimum eigenvalue of the random Gaussian matrix is lower bounded by a
positive value. This statement is true by Theorem 1.2 from Rudelson &
Vershynin (2008) since here each row of
$\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}$ has the distribution $\sim
s\kappa_{s}\mathcal{N}(0,\mathbb{I}_{s})$. Thus the random Gaussian matrix is
positive definite w.h.p. Therefore, we have statement (3) since
$\sigma_{\min}[AB]\geq\sigma_{\min}[A]\cdot\sigma_{\min}[B]$ where $A,B$ are
conformable Hermitian matrices. (4) holds subsequently from Theorem 1.2 from
Rudelson & Vershynin (2008) by considering the random matrix $n\cdot N$ with
entries from the distribution $s\kappa_{s}\mathcal{N}(0,\mathbb{I}_{s^{2}})$.
(5) is based on the bounds of the eigenvalues of a RIP matrix $K$ with RIC
$\zeta_{s+1}$ and order $s+1$.

######  Lemma 20.

Consider the term $M_{2}=\mathbb{I}-(\mathbb{I}+NA^{-1})^{-1}NA^{-1}$. We show
that $||M_{2}||_{2}$ is upper bounded by $\kappa_{M_{2}}$ where
$N=\frac{\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}}{n}$ and
$A=\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n}$ (as terms represented
in Algorithm 3). Here,
$\kappa_{M_{2}}=1-\frac{\kappa_{M_{1}}\varepsilon_{\eta}\sqrt{s}\kappa_{s}}{(1+\zeta_{s+1})n}$
and
$\kappa_{M_{1}}=\frac{||v||_{2}}{1+\frac{C_{1}s^{3/2}\kappa_{s}}{n(1-\zeta_{s+1})}}=\kappa_{M_{1}}||v||_{2}$.
Furthermore, $||M_{1}||=||(\mathbb{I}+NA^{-1})v||_{2}$ is upper bounded by
$\kappa_{M_{1}}$. Here, $\varepsilon_{\eta},C_{1}$ depend (polynomially) on
the sub-Gaussian moment bounds of $\mathcal{N}(0,1)$ (Rudelson & Vershynin,
2008).

Proof:

$\displaystyle||M_{2}||_{2}$
$\displaystyle=||\mathbb{I}-(\mathbb{I}+NA^{-1})^{-1}NA^{-1}||_{2}$
$\displaystyle\overset{(1)}{=}1-||(\mathbb{I}+NA^{-1})^{-1}NA^{-1}||_{2}$
$\displaystyle\overset{(2)}{\leq}1-\kappa_{M_{1}}||NA^{-1}||_{2}$
$\displaystyle\overset{(3)}{\leq}1-\frac{\kappa_{M_{1}}\varepsilon_{\eta}s\kappa_{s}}{n\sqrt{s}}||A^{-1}||_{2}$
$\displaystyle\overset{(4)}{\leq}1-\frac{\kappa_{M_{1}}\varepsilon_{\eta}\sqrt{s}\kappa_{s}}{n(1+\zeta_{s+1})}$

Here statement (1) follows by $||k\mathbb{I}-K||_{2}=k-||K||_{2}$ where
$\mathbb{I}$ is the identity matrix and $K$ is a conformable matrix. (2) is
based on identifying a lower bound for $M_{1}$ (below). (3) employs theorem
1.2 from Rudelson & Vershynin (2008) where the rows in $n\cdot N$ have
distribution of $s\kappa_{s}\mathcal{N}(0,\mathbb{I}_{s})$. (4) uses the bound
on the eigen-values of the RIP matrix $A$. To compute the lower bound for
$||N||$ we consider some conformable vector $v$,

$\displaystyle||M_{1}||_{2}$
$\displaystyle=||(\mathbb{I}+NA^{-1})^{-1}v||_{2}$
$\displaystyle\overset{(1)}{\geq}\sigma_{\min}((\mathbb{I}+NA^{-1})^{-1})||v||_{2}$
$\displaystyle\overset{(2)}{\geq}\frac{||v||_{2}}{\sigma_{\max}(\mathbb{I}+NA^{-1})}$
$\displaystyle\overset{(3)}{=}\frac{||v||_{2}}{1+\sigma_{\max}(NA^{-1})}$
$\displaystyle\overset{(4)}{\geq}\frac{||v||_{2}}{1+\sigma_{\max}(N)\sigma_{\max}(A^{-1}))}$
$\displaystyle\overset{(5)}{\geq}\frac{||v||_{2}}{1+\frac{C_{1}\sqrt{s}s\kappa_{s}}{(n)(1-\zeta_{s+1})}}$
$\displaystyle\overset{(6)}{=}\frac{||v||_{2}}{1+\frac{C_{1}s^{3/2}\kappa_{s}}{n(1-\zeta_{s+1})}}=\kappa_{M_{1}}||v||_{2}$

(1)-(4) follow from simple matrix manipulations, (5) is based on Theorem 2.4
from (Rudelson & Vershynin, 2008) for matrix $n\cdot N$ with values
distributed as $s\kappa_{s}\mathcal{N}(0,1)$ and (6) is based on bounds on
eigenvalues of matrix $A^{-1}$.

######  Lemma 21.

We show that the terms 2-norm of the terms
$L_{1}=\frac{X_{S_{\mathcal{A}}}^{T}(X_{S_{\mathcal{A}}^{c}}\alpha_{S_{*}}^{c}+\epsilon)}{n}$
and $L_{2}=\frac{X_{S_{\mathcal{A}}}^{T}y}{n}$ required in our analysis are
upper-bounded as expressed below. Essentially, $L_{1}$ considers the value of
$y$ after multiplying it by the projection matrix that removes certain columns
of $X_{S_{\mathcal{A}}}$ (due to orthogonality between these columns and the
projection matrix).

$\displaystyle||L_{1}||_{2}$
$\displaystyle\leq\zeta_{s+1}||\alpha_{S_{*}}^{c}||_{2}+\sqrt{1+\zeta_{s+1}}\kappa_{\epsilon}$
$\displaystyle||L_{2}||_{2}$
$\displaystyle\leq\sqrt{1+\zeta_{s+1}}\cdot\Big{(}||\alpha_{S_{*}}||_{2}+\sqrt{1+\zeta_{s+1}}\kappa_{\epsilon}\Big{)}$

Proof: Consider the analysis for $L_{1}$,

$\displaystyle||L_{1}||_{2}$
$\displaystyle=||\frac{X_{S_{\mathcal{A}}}^{T}(X_{S_{\mathcal{A}}^{c}}\alpha_{S_{*}}^{c}+\epsilon)}{n}||_{2}$
$\displaystyle\overset{(1)}{\leq}||\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}^{c}}\alpha_{S_{*}}^{c}}{n}||_{2}+||\frac{X_{S_{\mathcal{A}}}^{T}e}{n}||_{2}$
$\displaystyle\overset{(2)}{\leq}(\zeta_{s+1})||\alpha_{S_{*}}^{c}||_{2}+\sqrt{1+\zeta_{s+1}}\kappa_{\epsilon}$

(1) follows by the triangle inequality and the RIP definition for matrices
$X_{S_{\mathcal{A}}}^{T}/\sqrt{n}$, $X_{S_{*}}/\sqrt{n}$ respectively. By the
$L_{2}$-$L_{\infty}$ inequality we further have
$||x||_{2}\leq\sqrt{d}||x||_{\infty}$ where $x\in\mathbb{R}^{d}$. Clearly,
with Lemma 36 we then have $||e||_{2}\leq\sqrt{n}\kappa_{\epsilon}$ w.h.p..
where $\kappa_{\epsilon}=\sigma_{\epsilon}\sqrt{2\log 2/p_{0}}$ where
$e_{i}\sim\mathcal{N}(0,\sigma_{\epsilon}^{2})$ and $p_{0}$ is determined by
our probability requirements. Now, for $L_{2}$,

$\displaystyle||L_{2}||_{2}$
$\displaystyle=||\frac{X_{S_{\mathcal{A}}}^{T}y}{n}||_{2}$
$\displaystyle=||\frac{X_{S_{\mathcal{A}}}^{T}(X_{S_{*}}\alpha_{S_{*}}+\epsilon)}{n}||_{2}$
$\displaystyle\overset{(1)}{\leq}||\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{*}}\alpha_{S_{*}}}{n}+\frac{X_{S_{\mathcal{A}}}^{T}e}{n}||_{2}$
$\displaystyle\overset{(2)}{\leq}(1+\zeta_{s+1})||\alpha_{S_{*}}||_{2}+\sqrt{1+\zeta_{s+1}}\kappa_{\epsilon}$

(1) follows since $y=X_{S_{*}}\alpha_{S_{*}}+\epsilon$. (2) follows by RIP
definition and Lemma 18 for the first term and Lemmas 36,18 for the error
term.

###  D.1 Bounding the noise terms in SPriFed-OMP

######  Lemma 22.

[Bound small inverse term] We show that for a small value $0<t<1$,
$\frac{1}{1-t}\leq 1+\nu t$ if $\nu=\frac{1}{1-B_{t}}$ where $B_{t}$ is an
upper bound of $t$. Furthermore, $\nu\in(1,2)$.

Proof: Suppose, we assume that $\frac{1}{1-t}\leq 1+\nu t$ holds. Then, we can
write the inequality simply,

$\displaystyle 0\leq(\nu-1)t-\nu t^{2}$ $\displaystyle\rightarrow
t\leq\frac{\nu-1}{\nu}$

Setting, $B_{t}=\frac{\nu-1}{\nu}$, we obtain $\nu=\frac{1}{1-B_{t}}$. We
note, that for the second statement to hold, we require that $\nu-1>0$ and for
$\nu\geq 2$ the statement holds trivially.

####  D.1.1 Note on the difference between the basis columns and the non-basis
columns

We present the analysis for noise terms when only a single feature $j$ is
picked each time (Lemmas 23-27). However, we note that the same result holds
when we pick the set $S_{\mathcal{A}}^{c}$ instead of $j$ alone. This is true
since the total cardinality of the first matrix (either for the noise or
signal matrix) is less than $s+1$. Therefore, results from RIP and other
singular value bounds apply to this case as well. This is particularly
important for true basis columns. However, for the non-basis columns, this is
not useful as the total cardinality of the non-basis set and the SPriFed-OMP
selected feature set will easily exceed $s$ (closer to $p$).

######  Lemma 23.

The noise terms

$\displaystyle||\eta_{1}^{j}||$
$\displaystyle=||X^{T}_{j}X_{S_{\mathcal{A}}}(X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}})^{-1}$
$\displaystyle\cdot[\mathbb{I}{+}\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}(X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}})^{-1}]^{-1};j\in
S_{*}$ $\displaystyle||(\eta_{1}^{j})^{{}^{\prime}}||$
$\displaystyle=||X^{T}_{j}X_{S_{\mathcal{A}}}(X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}})^{-1}$
$\displaystyle\cdot[\mathbb{I}{+}\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}(X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}})^{-1}]^{-1};j\notin
S_{*}$

(Theorem 7) are upper-bounded such that,

$\displaystyle||(\eta_{1}^{j})^{{}^{\prime}}||_{2}$ $\displaystyle\leq
s^{2}\zeta_{s+1}\kappa_{s}\sqrt{2\log(\frac{s^{2}}{p_{b}})}(1+\nu\zeta_{s+1})^{2}\kappa_{M}$
$\displaystyle\cdot(\sqrt{1+\zeta_{s+1}}||\alpha_{S_{*}}||_{2}+\kappa_{\epsilon})$
$\displaystyle||\eta_{1}^{j}||_{2}$ $\displaystyle\leq
s^{2}\kappa_{s}\sqrt{2\log(\frac{s^{2}}{p_{b}})}(1+\nu\zeta_{s+1})^{3}\kappa_{M}$
$\displaystyle\cdot(\sqrt{1+\zeta_{s+1}}||\alpha_{S_{*}}||_{2}+\kappa_{\epsilon})$

Here, terms are borrowed from Theorem 7, Lemma 19) and $p_{b}$ is base failure
probability defined in Lemma 36.

Proof: We provide the proof for the $||(\eta_{1}^{j})^{{}^{\prime}}||_{2}$
term. However, proving the inequality for the $||\eta_{1}^{j}||_{2}$ is
straightforward.

$\displaystyle||(\eta_{1}^{j})^{{}^{\prime}}||_{2}$
$\displaystyle=||X^{T}_{j}X_{S_{\mathcal{A}}}(X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}})^{-1}[\mathbb{I}+\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}(X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}})^{-1}]^{-1}$
$\displaystyle\cdot\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}(X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}})^{-1}X_{S_{\mathcal{A}}}^{T}y||_{2}$
$\displaystyle=n||\frac{X^{T}_{S^{c}_{\mathcal{A}}}X_{S_{\mathcal{A}}}}{n}(\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n})^{-1}$
$\displaystyle\cdot[\mathbb{I}+\frac{\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}}{n}(\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n})^{-1}]^{-1}$
$\displaystyle\cdot\frac{\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}}{n}(\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n})^{-1}\frac{X_{S_{\mathcal{A}}}^{T}y}{n}||_{2}$
$\displaystyle\overset{(1)}{\leq}n\frac{\zeta_{s+1}}{(1-\zeta_{s+1})^{2}}||M||_{2}||\frac{\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}}{n}||_{2}L_{1}$
$\displaystyle\overset{(2)}{\leq}n\zeta_{s+1}||\frac{\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}}{n}||_{2}(1+\nu\zeta_{s+1})^{2}\kappa_{M}$
$\displaystyle\cdot(\sqrt{1+\zeta_{s+1}}||\alpha_{S_{*}}||_{2}+\kappa_{\epsilon})$
$\displaystyle\overset{(3)}{\leq}s^{2}\zeta_{s+1}\kappa_{s}\sqrt{2\log(\frac{s^{2}}{p_{b}})}(1+\nu\zeta_{s+1})^{2}\kappa_{M}$
$\displaystyle\cdot(\sqrt{1+\zeta_{s+1}}||\alpha_{S_{*}}||_{2}+\kappa_{\epsilon})$

Here, (1) follows from Lemma 17, for matrices $\frac{X^{T}_{j}}{\sqrt{n}}$ and
$\frac{X_{S_{\mathcal{A}}}}{\sqrt{n}}$ and bounded eigenvalues for
$||(\frac{X^{T}_{S_{\mathcal{A}}}X_{S_{\mathcal{A}}}}{n})^{-1}||_{2}$. The
middle and last terms are handled by Lemmas 20 and 21 respectively. In (2), we
replace the bounds for the terms $M,L_{1}$ using Lemmas 20 and 21
respectively. We replace the denominator by using the bound in Lemma 22.
Finally, (3) follows from simply noticing that each element in
$\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}$ has distribution
$\mathcal{N}(0,s^{2}\kappa_{s}^{2})$, $L_{2}-L_{\infty}$ inequality and lemma
36. The bound for the noise term holds with probability close to $1$ of the
form $p_{b}=1-k_{1}e^{-k_{2}}\approx 1$ where $k_{1},k_{2}$ are positive
constants.

######  Lemma 24.

The noise terms

$\displaystyle||\eta_{31}^{j}||_{2}$
$\displaystyle=n||\frac{X^{T}_{j}X_{S_{\mathcal{A}}}}{n}(\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n})^{-1}\Bigg{(}\mathbb{I}-[\mathbb{I}+\frac{\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}}{n}(\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n})^{-1}]^{-1}$
$\displaystyle\cdot\frac{\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}}{n}(\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n})^{-1}\Bigg{)}\frac{\eta_{\gamma_{\mathcal{S_{\mathcal{A}}}}}}{n}||_{2};j\in
S_{*}$ $\displaystyle||(\eta_{31}^{j})^{{}^{\prime}}||_{2}$
$\displaystyle=n||\frac{X^{T}_{j}X_{S_{\mathcal{A}}}}{n}(\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n})^{-1}\Bigg{(}\mathbb{I}-[\mathbb{I}+\frac{\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}}{n}$
$\displaystyle\cdot(\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n})^{-1}]^{-1}$
$\displaystyle\cdot\frac{\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}}{n}(\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n})^{-1}\Bigg{)}\frac{\eta_{\gamma_{\mathcal{S_{\mathcal{A}}}}}}{n}||_{2};j\notin
S_{*}$

is shown to be upper bounded such that,

$\displaystyle||(\eta_{31}^{j})^{{}^{\prime}}||_{2}$ $\displaystyle\leq
s\kappa_{s}\zeta_{s+1}\sqrt{2\log(\frac{2s}{p_{b}})}(1+\nu\zeta_{s+1})\kappa_{M_{2}}$
$\displaystyle||\eta_{31}^{j}||_{2}$ $\displaystyle\leq
s\kappa_{s}\sqrt{2\log(\frac{2s}{p_{b}})}(1+\nu\zeta_{s+1})^{2}\kappa_{M_{2}}$

The terms follow the notation in Algorithm 3, Theorem 7, Lemmas 22, 20.
$p_{b}$ is the base failure probability defined in Lemma 36.

Proof: We first demonstrate the proof for the case when $j\notin S_{*}$. The
proof for the $j\in S_{*}$ case follows simply by replacing the first
inequality involving the term $X_{j}^{T}X_{S_{\mathcal{A}}}$.

$\displaystyle||(\eta_{31}^{j})^{{}^{\prime}}||_{2}$
$\displaystyle=n||\frac{X^{T}_{j}X_{S_{\mathcal{A}}}}{n}(\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n})^{-1}\Bigg{(}\mathbb{I}-[\mathbb{I}+\frac{\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}}{n}$
$\displaystyle\cdot(\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n})^{-1}]^{-1}$
$\displaystyle\cdot\frac{\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}}{n}(\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n})^{-1}\Bigg{)}\frac{\eta_{\gamma_{\mathcal{S_{\mathcal{A}}}}}}{n}||_{2}$
$\displaystyle\overset{(1)}{\leq}\frac{n\zeta_{s+1}}{(1-\zeta_{s+1})}||M_{2}||_{2}||\frac{\eta_{\gamma_{\mathcal{S_{\mathcal{A}}}}}}{n}||_{2}$
$\displaystyle\overset{(2)}{\leq}\frac{n(\zeta_{s+1})\kappa_{M_{2}}}{(1-\zeta_{s+1})}||\frac{\eta_{\gamma_{\mathcal{S_{\mathcal{A}}}}}}{n}||_{2}$
$\displaystyle\overset{(3)}{\leq}\frac{s\kappa_{s}\sqrt{2\log(\frac{2s}{p_{b}})}\zeta_{s+1}\kappa_{M_{2}}}{(1-\zeta_{s+1})}$
$\displaystyle\overset{(4)}{\leq}s\kappa_{s}\zeta_{s+1}\sqrt{2\log(\frac{2s}{p_{b}})}(1+\nu\zeta_{s+1})\kappa_{M_{2}}$

Here, (1) follows from RIP properties in Lemmas 16-17 (since $j\notin
S_{\mathcal{A}})$ and $M_{2}$ is defined in Lemma 20. (2) follows from the
bound in Lemma 20. (3) follows from Lemma 36, $L_{2}$-$L_{\infty}$ inequality
and the Gaussian Mechanism for differential privacy. Here,
$(\eta_{\gamma_{S_{\mathcal{A}}}})^{j}\sim\mathcal{N}(0,\kappa_{s}^{2}s)$. And
$\kappa_{s}$ is defined in Algorithm 3. (4) is a direct application of Lemma
22.

######  Lemma 25.

The noise terms $||\eta_{32}^{j}||_{2},||(\eta_{32}^{j})^{{}^{\prime}}||_{2}$
which both have the form
$n||\frac{\eta_{\beta_{j,S_{\mathcal{A}}}}}{n}(\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n})^{-1}\Big{(}\mathbb{I}-[\mathbb{I}+\frac{\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}}{n}(\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n})^{-1}]^{-1}\frac{\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}}{n}(\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n})^{-1}\Big{)}\frac{X_{S_{\mathcal{A}}}^{T}y}{n}||_{2}$
for sets $j\in S_{*}$ and $j\notin S_{*}$ respectively are both upper bounded
by
$\sqrt{ps}\kappa_{p}\sqrt{2\log(\frac{2p}{p_{b}})}\kappa_{M_{2}}\Big{(}\sqrt{1+\zeta_{s+1}}||\alpha_{S_{*}}||_{2}+\kappa_{\epsilon}\Big{)}(1+\nu\zeta_{s+1}))$.
The notation is borrowed from Algorithm 3 and Theorem 7, Lemmas 22, 20.
$p_{b}$ is the base failure probability defined in Lemma 36.

$\displaystyle||\eta_{32}^{j}||_{2}$
$\displaystyle=n||\frac{\eta_{\beta_{j,S_{\mathcal{A}}}}}{n}(\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n})^{-1}\Big{(}\mathbb{I}-[\mathbb{I}+\frac{\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}}{n}$
$\displaystyle\cdot(\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n})^{-1}]^{-1}$
$\displaystyle\cdot\frac{\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}}{n}(\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n})^{-1}\Big{)}\frac{X_{S_{\mathcal{A}}}^{T}y}{n}||_{2}$
$\displaystyle\overset{(1)}{\leq}n||\frac{\eta_{\beta_{j,S_{\mathcal{A}}}}}{n}||_{2}\frac{||M_{2}||_{2}||L_{2}||_{2}}{1-\zeta_{s+1}}$
$\displaystyle\overset{(2)}{\leq}n||\frac{\eta_{\beta_{j,S_{\mathcal{A}}}}}{n}||_{2}\frac{\kappa_{M_{2}}\Big{(}\sqrt{1+\zeta_{s+1}}||\alpha_{S_{*}}||_{2}+\kappa_{\epsilon}\Big{)}}{1-\zeta_{s+1}}$
$\displaystyle\overset{(3)}{\leq}n||\frac{\eta_{\beta_{j,S_{\mathcal{A}}}}}{n}||_{2}\kappa_{M_{2}}\Big{(}\sqrt{1+\zeta_{s+1}}||\alpha_{S_{*}}||_{2}+\kappa_{\epsilon}\Big{)}$
$\displaystyle\cdot(1+\nu\zeta_{s+1}))$
$\displaystyle\overset{(4)}{\leq}n\sqrt{s}||\frac{\eta_{\beta_{j,S_{\mathcal{A}}}}}{n}||_{\infty}\kappa_{M_{2}}\Big{(}\sqrt{1+\zeta_{s+1}}||\alpha_{S_{*}}||_{2}+\kappa_{\epsilon}\Big{)}$
$\displaystyle\cdot(1+\nu\zeta_{s+1}))$
$\displaystyle\overset{(5)}{\leq}\sqrt{ps}\kappa_{p}\sqrt{2\log(\frac{2p}{p_{b}})}\kappa_{M_{2}}\Big{(}\sqrt{1+\zeta_{s+1}}||\alpha_{S_{*}}||_{2}+\kappa_{\epsilon}\Big{)}$
$\displaystyle\cdot(1+\nu\zeta_{s+1}))$

Statement (1) follows by substituting for the middle and last term (from
Lemmas 20, 21 resp.), separating the noise vector and bounding the singular
values of $\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n})^{-1}$. (2)
follows by bounds in Lemmas 20, 21. (3) is a direct application of inverting
the denominator term with Lemma 22. (4) uses the $L_{2}-L_{\infty}$
inequality. The final statement uses Lemma 36, distribution of the noise
matrix and bounds it for the appropriate set and failure probability as shown
in Lemma 36.

######  Lemma 26.

The noise terms $||\eta_{33}^{j}||_{2},||(\eta_{33}^{j})^{{}^{\prime}}||_{2}$
over sets $j\in S_{*}$ and $j\notin S_{*}$ respectively both have the form

$\displaystyle
n||\frac{\eta_{\beta_{j,S_{\mathcal{A}}}}}{n}(\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n})^{-1}\Big{(}\mathbb{I}-[\mathbb{I}+\frac{\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}}{n}(\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n})^{-1}]^{-1}$
$\displaystyle\cdot\frac{\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}}{n}(\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n})^{-1}\Big{)}\frac{\eta_{\gamma_{\mathcal{S_{\mathcal{A}}}}}}{n}||_{2}$

are both upper bounded by
$2\kappa_{M_{2}}(1+\nu\zeta_{s+1})\frac{s\kappa_{s}\kappa_{p}}{\kappa_{n}}\sqrt{\log(\frac{2p}{p_{b}})\log(\frac{2s}{p_{b}})}$.
The notation follows from Algorithm 3 and Theorem 7, Lemmas 22, 20. $p_{b}$ is
the base failure probability defined in Lemma 36.

$\displaystyle||\eta_{33}^{j}||_{2}$
$\displaystyle=n||\frac{\eta_{\beta_{j,S_{\mathcal{A}}}}}{n}(\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n})^{-1}\Big{(}\mathbb{I}-[\mathbb{I}+\frac{\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}}{n}(\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n})^{-1}]^{-1}\Big{)}$
$\displaystyle\cdot\frac{\eta_{\beta_{S_{\mathcal{A}},S_{\mathcal{A}}}}}{n}(\frac{X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}}{n})^{-1}\Big{)}\frac{\eta_{\gamma_{\mathcal{S_{\mathcal{A}}}}}}{n}||_{2}$
$\displaystyle\overset{(1)}{\leq}n\frac{||M_{2}||_{2}}{1-\zeta_{s+1}}||\frac{\eta_{\beta_{j,S_{\mathcal{A}}}}}{n}||_{2}||\frac{\eta_{\gamma_{\mathcal{S_{\mathcal{A}}}}}}{n}||_{2}$
$\displaystyle\overset{(2)}{\leq}n\kappa_{M_{2}}(1+\nu\zeta_{s+1})\sqrt{s}||\frac{\eta_{\beta_{j,S_{\mathcal{A}}}}}{n}||_{\infty}\sqrt{s}||\frac{\eta_{\gamma_{\mathcal{S_{\mathcal{A}}}}}}{n}||_{\infty}$
$\displaystyle\overset{(3)}{\leq}\kappa_{M_{2}}(1+\nu\zeta_{s+1})\frac{\sqrt{ps}s}{n}\kappa_{s}\kappa_{p}2\sqrt{\log(\frac{2p}{p_{b}})\log(\frac{2s}{p_{b}})}$
$\displaystyle\overset{(4)}{=}2\kappa_{M_{2}}(1+\nu\zeta_{s+1})\frac{s\kappa_{s}\kappa_{p}}{\kappa_{n}}\sqrt{\log(\frac{2p}{p_{b}})\log(\frac{2s}{p_{b}})}$

Statement (1) follows from bounding the eigenvalues of
$X_{S_{\mathcal{A}}}^{T}X_{S_{\mathcal{A}}}{n})^{-1}$, substituting the value
of the middle term as $M_{2}$ (from Lemma 20) and separating the noise
2-norms. (2) leverages the $L_{2}-L_{\infty}$ inequality for the noise
2-norms, bounds $M_{2}$ with Lemma 20 and inverts the denominator with the
simple Lemma 22. (3) follows from noting the variances of the noise matrices
as being proportional to $p$ and $s$ and using the maximum concentration bound
(Lemma 36) to pick an appropriate value. $\kappa_{s}$ and $\kappa_{p}$ are
picked by the privacy mechanism in Algorithm 3. Depending on the set we are
maximizing over for $j$ we can pick $d=s$ (basis columns) or $d=p-s$ (non-
basis columns). $p_{b}$ is some appropriate base probability that is generally
quite small. (4) follows by substituting $n=\sqrt{ps}\kappa_{n}$.

######  Lemma 27.

The noise term
$||\frac{\eta^{j}_{2}}{n}||_{2}=||\frac{(\eta_{\gamma_{0}^{-}})_{j}}{n}||_{2}$
is upper bounded by
$\frac{\sqrt{2\log(\frac{2p}{p_{b}})}\kappa_{p}}{\sqrt{s}\kappa_{n}}$. Given,
$j\in T,d=|T|$ and $p_{b}$ is some base probability that is quite small.

Proof:

$\displaystyle||\frac{\eta^{j}_{2}}{n}||_{2}$
$\displaystyle=||\frac{(\eta_{\gamma_{0}^{-}})_{j}}{n}||_{2}\leq\frac{\sqrt{p}}{n}\kappa_{p}\sqrt{2\log(\frac{2p}{p_{b}})}=\frac{\sqrt{2\log(\frac{2p}{p_{b}})}\kappa_{p}}{\sqrt{s}\kappa_{n}}$

where the last inequality holds because of Lemma 36 and since
$(\eta_{\gamma_{0}^{-}})_{j}\sim\mathcal{N}(0,\kappa_{p}^{2}p)$ where
$\kappa_{p}$ is picked by the privacy mechanism in Algorithm 3.

##  Appendix E Proof of Main Theorem 7

We will use induction to prove our main result as noted in section 6.

Step 1: We show that we pick the correct basis in step 1 with high probability
(w.h.p.). Here, the correct basis implies that we choose the column index
belonging to $\mathcal{S_{*}}$. First step from Algorithm 3 is written as
$j\leftarrow\text{argmax}_{j\in\Omega}|(\gamma_{0}^{-})_{j}|$. Our strategy is
to find the maximum correlation’s lower and upper bound under the assumptions
$j\in S_{*}$ and $j\notin S_{*}$, respectively. Denote these correlations $C$
and $C^{{}^{\prime}}$ respectively and let their lower bound and upper bound
be $C_{*}$ and $C^{{}^{\prime}}_{*}$ respectively. If
$C_{*}^{{}^{\prime}}<C_{*}$ then the correct basis (from $S_{*}$) is chosen at
each step. This statement is true since $C_{*}^{{}^{\prime}}<C_{*}$ implies
that the upper bound of the maximum correlation of non-support columns is
smaller than the lower bound of the maximum correlation support columns. Thus,
we first attempt to find the values $C_{*},C^{{}^{\prime}}_{*}$ so as to
equate the inequality above.

Thus, the lower bound of the maximum correlation of support columns is given
by,

$\displaystyle C$ $\displaystyle=\max_{j\in S_{*}}|(\gamma_{0}^{-})_{j}|$
$\displaystyle=\max_{j\in S_{*}}|X_{j}^{T}y+\eta_{(\gamma_{0}^{-})_{j}}|$
$\displaystyle\geq\frac{||X_{S_{*}}^{T}y+\eta_{(\gamma_{0}^{-})_{j}}||_{2}}{\sqrt{s}}\text{
(by $L_{2},L_{\infty}$ ineq.)}$
$\displaystyle\geq\frac{||X_{S_{*}}^{T}y||_{2}}{\sqrt{s}}-\max_{j\in
S_{*}}|\eta_{(\gamma_{0}^{-})_{j}}|\text{ (by reverse triangle inequality)}$

Suppose events $E_{1},E_{2},...,E_{s}$ each occur w.p. $(1-p_{a})$. Event
$E_{k},k\in[s]$ indicates that $|\eta_{(\gamma_{0}^{-})_{k}}|\leq B_{a}$ for
some $B_{a}>0$. This could be rewritten as
$\text{Pr}[E_{k}]=\text{Pr}[|\eta_{(\gamma_{0}^{-})_{k}}|\leq B_{a}]=1-p_{a}$.
Invoking Lemma 32 and since
$|\eta_{(\gamma_{0}^{-})_{k}}|\sim\mathcal{N}(0,\sigma^{2}_{0})$ we can write
$p_{a}\triangleq 2e^{-\frac{B_{a}}{2\sigma_{0}^{2}}}\rightarrow
B_{a}=\sigma_{0}\sqrt{2\text{log}(\frac{2}{p_{a}})}=\sqrt{2\text{log}(\frac{2}{p_{a}})}\frac{\kappa_{p}}{\epsilon_{step}}\sqrt{p}=\kappa_{a}\cdot\sqrt{p}$
where
$\kappa_{a}=\sqrt{2\text{log}(\frac{2}{p_{a}})}\frac{\kappa_{p}}{\epsilon_{step}}$.
To bound the maximum of the noise terms, we require that all events
$E_{k},k\in[s]$ hold simultaneously. By the union bound (Lemma 33),
$\text{Pr}[\max_{k\in
S_{*}}|\eta_{(\gamma_{0}^{-})_{k}}|>B_{a}]=\text{Pr}[\cup_{k\in
S_{*}}E_{k}^{c}]\leq\sum\limits_{k\in S_{*}}\text{Pr}[E_{k}^{c}]=s\cdot
p_{a}$. Therefore, $\text{Pr}[\max_{k\in
S_{*}}|\eta_{(\gamma_{0}^{-})_{k}}|\leq B_{a}]\geq 1-s\cdot p_{a}$. And so
w.p. $1-s\cdot p_{a}$ we can say that,

$\displaystyle C$
$\displaystyle\geq\frac{||X_{S_{*}}^{T}y||_{2}}{\sqrt{s}}-\kappa_{a}\sqrt{p}$
$\displaystyle=\frac{||X_{S_{*}}^{T}X_{S_{*}}\alpha_{S_{*}}+X_{S_{*}}^{T}e||_{2}}{\sqrt{s}}-\kappa_{a}\sqrt{p}$