# Private and polynomial time algorithms for learning Gaussians and beyond

Hassan Ashtiani McMaster University<EMAIL_ADDRESS>Hassan Ashtiani is
also affiliated with Vector institute and was supported by an NSERC Discovery
Grant. Christopher Liaw Google<EMAIL_ADDRESS>Part of this work was done
while at the University of Toronto.

###### Abstract

We present a fairly general framework for reducing $(\varepsilon,\delta)$
differentially private (DP) statistical estimation to its non-private
counterpart. As the main application of this framework, we give a polynomial
time and $(\varepsilon,\delta)$-DP algorithm for learning (unrestricted)
Gaussian distributions in $\mathbb{R}^{d}$. The sample complexity of our
approach for learning the Gaussian up to total variation distance $\alpha$ is
$\widetilde{O}(d^{2}/\alpha^{2}+d^{2}\sqrt{\ln(1/\delta)}/\alpha\varepsilon+d\ln(1/\delta)/\alpha\varepsilon)$
matching (up to logarithmic factors) the best known information-theoretic
(non-efficient) sample complexity upper bound due to Aden-Ali, Ashtiani, and
Kamath [1]. In an independent work, Kamath, Mouzakis, Singhal, Steinke, and
Ullman [25] proved a similar result using a different approach and with
$O(d^{5/2})$ sample complexity dependence on $d$.

As another application of our framework, we provide the first polynomial time
$(\varepsilon,\delta)$-DP algorithm for robust learning of (unrestricted)
Gaussians with sample complexity $\widetilde{O}(d^{3.5})$. In another
independent work, Kothari, Manurangsi, and Velingker [30] also provided a
polynomial time $(\varepsilon,\delta)$-DP algorithm for robust learning of
Gaussians with sample complexity $\widetilde{O}(d^{8})$.

Accepted for presentation at the Conference on Learning Theory (COLT) 2022
[4].

##  1 Introduction

inline,color=red!30]Chris: Cite Pravesh Kothari paper. Learning a multivariate
Gaussian distribution with respect to the total variation (TV) distance is one
of the most basic tasks in statistical estimation. It is folklore that simply
taking the empirical mean and covariance matrix achieves the optimal sample
complexity of $\Theta(d^{2}/\alpha^{2})$ for learning _unbounded_ Gaussians in
$\mathbb{R}^{d}$ up to TV distance $\alpha$ (see, for example, [3, Theorem
C.1]). However, this simple approach can compromise the privacy of the
individuals/entities whose data has been used in the estimation procedure.
Differential privacy [19] offers a formal framework to study and guarantee
privacy in such estimation tasks. The central question that we aim to address
is whether unbounded Gaussians can be learned _privately_ and _efficiently_.

For the strong notion of _pure_ ($\varepsilon,0$)-differential privacy,
learning unbounded Gaussians is impossible. In fact, even learning the mean of
univariate Gaussians with unit variances requires some boundedness assumptions
on the parameters of the Gaussian [29, Theorem 1.4]. This motivates the study
of learning Gaussians under the _approximate_
($\varepsilon,\delta$)-differential privacy model [17] (or approximate-DP for
short).

Previous work has established the polynomial-time learnability of univariate
Gaussians (with unbounded mean and variance) in the approximate-DP model. The
seminal work of Karwa and Vadhan [29] uses stability-based histograms [8] to
find a bound on the parameters; once the problem is reduced to the bounded
case, the Laplace mechanism [19] offers a private way to release a noisy—but
sufficiently accurate—estimate of the parameters of the Gaussian111Although
the method of [29] works with a histogram with an infinite number of bins, it
can be easily implemented efficiently by only constructing the bins that are
non-empty (i.e., lazy initialization)..

The problem is more difficult in the multivariate setting. Although one could
apply the method of [29] along each axis, this would only yield an accurate
solution for axis-aligned Gaussians, i.e. Gaussians whose covariance matrix is
a diagonal matrix.

The first result for multivariate Gaussians is due to [24]. They gave an
algorithm for learning multivariate Gaussians in the concentrated DP model222
Concentrated DP is a weaker notion of privacy than pure DP but stronger than
approx DP [21, 9]. . However, their algorithm requires prior knowledge of the
condition number of the covariance matrix and its sample complexity depends
polylogarithmically on the condition number. In a follow-up work, [5] proposed
the CoinPress algorithm which is a simple and practical algorithm for learning
multivariate Gaussians in the concentrated DP model with more or less the same
sample complexity.

More recently, Aden-Ali, Ashtiani, and Kamath [1] established a sharp upper
bound on the sample complexity of learning Gaussians in the approximate-DP
model, demonstrating that the dependence on the condition number (or the size
of the parameters) can be totally avoided. Unfortunately, the approach of [1]
is non-constructive and does not yield a computationally efficient method. In
fact, it only shows the _existence_ of a “locally small” cover for Gaussians
in TV distance, which yields an upper bound on the sample complexity through
the use of the private hypothesis selection method [7]. This raises the
following concrete question.

> Is there an $(\varepsilon,\delta)$-DP and polynomial time333 By polynomial-
> time, we mean $\operatorname{poly}(d,1/\varepsilon,1/\alpha,\log(1/\delta))$
> where $\alpha$ is the accuracy parameter in TV distance.  algorithm for
> learning unbounded Gaussian distributions over $\mathbb{R}^{d}$?

More ambitiously, we can ask if there is an efficient algorithm for learning
high-dimensional Gaussians _robustly_ , where a small fraction of the samples
are arbitrarily corrupted by an adversary. The breakthrough of [12, 32]
answered this question positively. In particular, the work of [12] offers a
polynomial time algorithm for robustly learning Gaussians with respect to the
TV distance. This raises the following natural question.

> Is there an $(\varepsilon,\delta)$-DP, _robust_ , and polynomial time
> algorithm for learning unbounded Gaussian distributions over
> $\mathbb{R}^{d}$?

###  1.1 Main Results

In this paper, we devise a simple and general algorithmic technique for
private statistical estimation. Our first application is an efficient
algorithm for learning a Gaussian distribution in $\mathbb{R}^{d}$.

######  Theorem 1.1 (Informal).

There is an $(\varepsilon,\delta)$-DP algorithm that takes
$m=\widetilde{O}\left(\frac{d^{2}}{\alpha^{2}}+\frac{d^{2}\sqrt{\ln{1/\delta}}}{\alpha\varepsilon}+\frac{d\ln(1/\delta)}{\alpha\varepsilon}\right)$
i.i.d. samples from an arbitrary Gaussian distributions over $\mathbb{R}^{d}$,
runs in time polynomial in $m$, and with high probability outputs a Gaussian
whose TV distance to the true distribution is smaller than $\alpha$.

The formal version of the theorem is stated in Theorem 5.7. The sample
complexity of this approach matches the best (non-efficient) learning method
[1] up to logarithmic factors. Moreover, [24, Theorem 56] and [29, Theorem
1.4] show sample complexity lower bounds of $\Omega(d/\alpha\varepsilon)$ and
$\Omega(\log(1/\delta)/\varepsilon)$ for private learning of spherical and
univariate Gaussians respectively. Finally, a lower bound of
$\tilde{\Omega}(d^{2}/\alpha^{2})$ is known for the non-private setting [3,
Theorem 1.2].

The main ingredient of proving Theorem 1.1 is a polynomial-time algorithm to
compute a preconditioner which brings the condition number of the covariance
matrix down to a constant. Once we have a preconditioner, we can then (mostly)
use known results to prove Theorem 1.1. The algorithm to compute a
preconditioner is extremely simple: it computes several candidate covariance
matrices from disjoint samples, computes a weighted average of the “good”
candidates, and then adds noise to resulting covariance matrix. Here, a
candidate covariance matrix is “good” if it is close to a large majority of
the other candidate covariance matrices; intuitively, this should also
indicate that it is close to the true covariance matrix. The first two steps
come directly from the framework discussed in Section 3 and the last step
(adding the noise) is discussed in Section 5.

In order to handle Gaussians with degenerate covariance matrices, we need a
private and efficient method for learning (exactly) the subspace of a
Gaussian. The recent work of Singhal and Steinke [39] offers a differentially
private but inefficient method for this problem. We develop an efficient
method for this problem which may be of independent interest (see Theorem 4.1
for a formal statement). We note that our result is somewhat incomparable to
that of Singhal and Steinke [39] since the assumptions are different.

######  Theorem 1.2 (Informal).

There is an $(\varepsilon,\delta)$-DP algorithm that takes
$m=O(d\ln(1/\delta)/\varepsilon)$ i.i.d. samples from a Gaussian distribution
in $\mathbb{R}^{d}$ and with probability 1, outputs the subspace that the
Gaussian is supported on. Furthermore, the running time of the algorithm is
$O(md)$.

In this paper, we prove Theorem 1.2 as an example of our technique from
Section 3. However, we remark that this is essentially equivalent to other
standard frameworks such as propose-test-release.

Next, we show how to learn Gaussians privately, robustly, and efficiently. An
$\alpha$-corrupted sample is a sample that is corrupted by an adversary as
follows. First, we draw $m$ i.i.d. samples from a Gaussian. Then, the
adversary chooses at most $\alpha m$ samples and modifies them arbitrarily. In
the robust case, we prove the following theorem; more details can be found in
7.5.

######  Theorem 1.3 (Informal).

There is a robust and $(\varepsilon,\delta)$-DP algorithm which given an
$\alpha$-corrupted sample of size
$m=\widetilde{\Omega}\left(d^{3.5}\ln(1/\delta)/\varepsilon\alpha^{3}\right)$
from $\mathcal{N}(\mu,\Sigma)$ with unknown $\mu\in\mathbb{R}^{d}$ and non-
singular $\Sigma\in\mathbb{R}^{d\times d}$, outputs $\widehat{\mu}$ and
$\widehat{\Sigma}$ such that the total variation distance between
$\mathcal{N}(\mu,\Sigma)$ and $\mathcal{N}(\widehat{\mu},\widehat{\Sigma})$ is
at $O(\alpha\ln(1/\alpha))$ with high probability. Furthermore, the running
time of the algorithm is $\operatorname{poly}(m)$.

Important ingredients of the above theorem are two subroutines for learning
the mean and the covariance robustly. The following result for mean estimation
is proved naturally in our framework, but we note [35, Theorem 7] has an
improved sample complexity of
$\widetilde{\Omega}(d/\alpha^{2}+d^{3/2}\ln(1/\delta)/\varepsilon\alpha)$ for
the same problem. Our paper is the first to solve the more difficult task of
private, robust, and efficient covariance estimation. The discussion of robust
learning of Gaussians can be found in Section 6 and Section 7.

######  Theorem 1.4 (Informal).

There is a robust and $(\varepsilon,\delta)$-DP algorithm which given an
$\alpha$-corrupted sample of size
$m=\widetilde{\Omega}\left(\frac{d^{3/2}\ln(1/\delta)}{\varepsilon\alpha^{2}}\right)$
from $\mathcal{N}(\mu,\Sigma)$ with unknown $\mu\in\mathbb{R}^{d}$ and
$0.5I_{d}\preceq\Sigma\preceq 2I_{d}$, outputs $\widetilde{\mu}$ such that
$\|\mu-\widetilde{\mu}\|_{2}=O(\alpha\sqrt{\log(1/\alpha)})$ with high
probability. Furthermore, the running time of the algorithm is
$\operatorname{poly}(m)$.

inline,color=green!30]Hassan: update the above sentence to reflect their new
result inline,color=red!30]Chris: Do we still want the robust mean estimation
or should we use Sewoong’s result? Their sample complexity is slightly better
than ours ($d/\alpha^{2}+d^{3/2}\ln(1/\delta)/\varepsilon\alpha$).

###  1.2 Techniques

inline,color=red!30]Chris: Please read.

We present a simple, efficient, and fairly general framework for private
statistical estimation which is inspired by the Propose-Test-Release framework
[18] and the Subsample-And-Aggregate framework [36]. At a high-level, our
framework works as follows. Suppose that we had a _non-private_ estimator
$\mathcal{A}$. Given a set of samples, we split the set into disjoint subsets
and run $\mathcal{A}$ on each of the subsets (in the spirit of Subsample-And-
Aggregate). We then privately check whether most of the (non-private) outcomes
are “close” to each other (in the spirit of Propose-Test-Release). If so, we
then aggregate the “good” results by compute a (clever) weighted average of
the “good” results and then releasing an estimate by noising the computed
average. The key is that the “good” results ought to be quite structured (for
example, coming from a small ball) and averaging helps by providing some
additional stability. This, in turns allows us to apply only a small amount of
noise. We remark that the aggregation component of our framework can be seen
as a special case of the FriendlyCore framework [41] which was developed
independently to our work (in particular, see Algorithm 4.1 in [41]). One
minor difference is that their aggregation is randomized while our aggregation
uses a weighted average.

In order for the above framework to work, we need two key properties. The
first is that if we run $\mathcal{A}$ on two disjoint i.i.d. samples and the
outputs of both are close to the true answer (for example, the true mean or
true covariance matrix) then the outputs are close to each other. This would
allow us to define “good” as simply outputs of $\mathcal{A}$ which are close
to other outputs of $\mathcal{A}$. The second property we require is that if
two weighted averages are computed using nearly the same weights then the
weighted average are close to each other. This allows us to assert that the
output is stable under small perturbations of the input, and, thus allows us
to add a small amount of noise to the weighted average. Both of these
properties are satisfied for normed spaces (such as the $\ell_{2}$ space).
Unfortunately, the space that we require for covariance matrices is _not_ a
normed space (it does not even satisfy the triangle inequality). However, it
does satisfy an _approximate_ triangle inequality for points that are
sufficiently close together. In Definition 3.1, we define the notion of a
_convex semimetric_ space, which, intuitively, slightly relaxes the
requirement of being a norm so that it is general enough to handle other
spaces such as the space of covariance matrices (with an appropriate “norm”).

The final part of the reduction is to output a “noisy” version of the average.
Note that at this point, we are essentially guaranteed that the average is
quite stable. In particular, changing a single input point results in very
small changes in the average. In Definition 3.3, we define the notion of a
_masking mechanism_. At a high-level, a masking mechanism $\mathcal{B}$ is a
mechanism such that $\mathcal{B}(Y_{1}),\mathcal{B}(Y_{2})$ are
indistinguishable provided $Y_{1},Y_{2}$ are sufficiently close. Note that a
masking mechanism, in and of itself, is not a differentially private algorithm
since it does not necessarily ensure that its inputs are close to each other.
For example, the Laplace and Gaussian mechanisms are masking mechanisms which
are often used to design differentially private algorithms. A key technical
ingredient for learning covariance matrices is to design an efficient masking
mechanism for the space of covariance matrices. We do this by viewing the
input covariance matrix as a vector of length $d^{2}$ and applying the
Gaussian mechanism but where the noise used for the Gaussian mechanism is
scaled by the input covariance matrix itself. A similar idea of using an
empirically rescaled Gaussian mechanism was also utilized by [6] to design an
algorithm for learning the mean under Mahalanobis distance without estimating
the covariance matrix.

###  1.3 Discussion of a concurrent result

In an independent work, Kamath, Mouzakis, Singhal, Steinke, and Ullman [25]
propose a different polynomial time algorithm for privately learning unbounded
high-dimensional Gaussian distributions. At a high-level, they devise an
iterative preconditioning method for the covariance matrix (in the spirit of
[24, 5]). To resolve the challenge of handling covariance matrices with
unbounded condition numbers, they propose a subspace learning method inspired
by the approach of [39]. The iterative approach to bound the ratio of
consecutive eigenvalues incurs additional cost of DP composition, resulting in
a term in the sample complexity of order $O(d^{5/2}/\varepsilon)$ rather than
$O(d^{2}/\varepsilon)$.

In contrast, our approach for preconditioning the covariance matrix is not
iterative, and in its core relies on a masking/noising mechanism for
covariance matrices. In particular, we are able to compute a preconditioner
with roughly $\widetilde{O}(d^{2}/\varepsilon)$ samples. Non-privately, $O(d)$
samples are needed to compute a preconditioner. It is an interesting open
problem to see if $O(d)$ samples suffices for the private setting as well. Our
general reduction also allows us to apply the result in other settings such as
_robust_ learning of Gaussians in polynomial time.

In another independent work, Kothari, Manurangsi, Velingker [30] proposed an
efficient algorithm for robust and differentially private estimation of the
mean, covariance matrix, and higher moments of distributions that satisfy
either a certifiable subgaussianity or certifiable hypercontractivity. Their
algorithm is designed using the sum-of-squares paradigm. In the case of
Gaussian distributions, their work obtains an efficient robust and
$(\varepsilon,\delta)$-DP algorithm for learning a Gaussian in
$\mathbb{R}^{d}$ with sample complexity
$\widetilde{\Omega}(d^{8}\ln^{4}(1/\delta)/\alpha^{4}\varepsilon^{4})$ where
$\alpha$ is the corruption parameter.

In yet another independent work, Tsfadia, Cohen, Kaplan, Mansour, and Stemmer
[41] introduced the FriendlyCore framework which allows for differentially
private aggregation. The aggregation component of our algorithm is quite
similar to that of FriendlyCore (see [41, Algorithm 4.3]) with a minor
difference being that we compute a weighted average instead of computing a
random “core”. Indeed, after the first version of our paper was posted on the
arXiv, Tsfadia, Cohen, Kaplan, Mansour, and Stemmer [41] showed how their
framework can also be used to learn the covariance matrix using the tools
developed in this paper (in particular, the masking mechanism for covariance
matrices in 5.1). For more details, see Appendix A of [41].

###  1.4 More on related work

inline,color=green!30]Hassan: add a citation to [diakonikolas2017being] and
maybe to [cheng2019high]. We don’t have a citation to [15] or [12, 32] in this
section either

inline,color=green!30]Hassan: add citations here and/or in the reduction to
FriendlyCore [41] and perhaps to the older GoodCenter [37]. Regarding the
aggregate approaches, perhaps we should cite this one too [40]? Most relevant
to our work are the results on learning Gaussians with respect to the TV
distance under $(\varepsilon,\delta)$-DP. Karwa and Vadhan [29] established a
polynomial-time and sample-efficient method for learning unbounded univariate
Gaussians under $(\varepsilon,\delta)$-DP. One can use this method for
learning multivariate Gaussians as well, by simply learning the one
dimensional projection of the Gaussians along each axis (and invoking the
composition for DP). While this approach works well for axis-aligned
Gaussians, for general Gaussians its sample complexity depends on the
condition number of the covariance matrix. Kamath, Li, Singhal, and Ullman
[24] consider the problem of privately learning high-dimensional
distributions. They provide a polynomial time algorithm for learning high-
dimensional Gaussians under the (stronger notion of) concentrated DP [21, 9].
The sample complexity of this approach improves over the multivariate version
of [29]; however, it still depends (logarithmically) on the condition number
of covariance matrix. On a follow up work, Biswas, Dong, Kamath, and Ullman
[5] propose the more practical method of CoinPress whose sample complexity
matches that of [24] for learning multivariate Gaussians. CoinPress is shown
to empirically work better than its counterparts—namely the multivariate
versions of [29, 16]—in high dimensions.

The best known lower bounds for learning $d$-dimensional Gaussians up to TV
distance $\alpha$ under $(\varepsilon,\delta)$-DP is
$\tilde{\Omega}(d/\alpha\varepsilon)$ [24]. In contrast, the sample complexity
of all the methods described above depends on the condition number of the
covariance matrix. Utilizing the private hypothesis selection framework
proposed in [7], Aden-Ali, Ashtiani, and Kamath [1] proved the first sample
complexity bound for privately learning general Gaussians that does not depend
on the condition number or the size of the parameters. The proof is, however,
non-constructive, as it invokes Zorn’s lemma to prove the existence of a
locally small cover for covariance matrices.

In the non-private setting, an $\alpha$-accurate estimation of the mean of a
$d$-dimensional Gaussian distribution with respect to the Mahalanobis distance
can be found using $O(d/\alpha^{2})$ samples. Achieving linear dependence on
$d$, however, is not easy in the private setting. Brown, Gaboardi, Smith,
Ullman, and Zakynthinou [6] propose a method for solving this problem under
$(\varepsilon,\delta)$-DP with linear dependence on $d$. Their method is,
however, not computationally efficient. One interesting aspect of [6] is the
use of a variant of the Gaussian mechanism that is scaled with the empirical
covariance. We also make use of a similar idea to noise the covariance matrix.

A related line of work is the problem of private mean estimation with respect
to the Euclidean distance. It is noteworthy that these results do not
translate to learning Gaussians with respect to the TV distance. Kamath,
Singhal, and Ullman [27] prove tight minimax rates for mean estimation for
distributions with bounded $k$-th moment. Liu, Kong, Kakade, and Oh [35]
proposed a robust, efficient, and private mean estimation method for sub-
Gaussian distributions with bounded mean.

Other related work include “universal” private estimators for mean, variance,
and scale under pure differential privacy [14], private Principal Component
Analysis [10, 46, 13, 34] and private learning of Gaussian mixture models [2,
26]. We refer the readers to [20, 43] for an introduction to differential
privacy, and to [28] for a survey on private statistical learning.

At a very high-level, the algorithms in this paper all follow the same
approach of privately aggregating the outputs of non-private algorithms. Our
approach can be viewed as combining the Subsample-And-Aggregate paradigm with
an aggregation scheme that is is similar to the FriendlyCore framework of
[41]. Older works have also made use of this blueprint including [37] in the
context of clustering and [40] in the context of statistical estimation.

##  2 Preliminaries

For a vector $v\in\mathbb{R}^{d}$, we use $\|v\|_{2}$ to denote the Euclidean
norm of the vector $v$. For a matrix $A\in\mathbb{R}^{d\times d}$, we use
$\|A\|_{F}$ to denote the Frobenius norm of $A$ and $\|A\|$ to denote the
operator norm of $A$. A matrix $A$ is positive semi-definite (resp. positive
definite), denoted $A\succeq 0$ (resp. $A\succ 0$) if $A$ is symmetric and all
its eigenvalues are non-negative (resp. positive). If $A\succeq 0$ then
$A^{1/2}$ denotes the unique positive semi-definite matrix satisfying
$A=(A^{1/2})^{2}$. Further, if $A\succ 0$ then $A^{-1/2}=(A^{1/2})^{-1}$. In
this paper, we also write $\mathcal{S}^{d}$ to denote the positive-definite
cone in $\mathbb{R}^{d\times d}$.

###  2.1 Probability Facts

In this paper, we use
$d_{\mathrm{TV}}\left({\mathcal{D}_{1}},{\mathcal{D}_{2}}\right)$ and
$D_{\mathrm{KL}}\left({\mathcal{D}_{1}}\parallel{\mathcal{D}_{2}}\right)$ to
denote the total variation distance and KL divergence between two
distributions $\mathcal{D}_{1},\mathcal{D}_{2}$.

######  Fact 2.1.

Let $\mu_{1},\mu_{2}\in\mathbb{R}^{d}$ and $\Sigma_{1},\Sigma_{2}\succ 0$.
Then

$\displaystyle
D_{\mathrm{KL}}\left({\mathcal{N}(\mu_{1},\Sigma_{1})}\parallel{\mathcal{N}(\mu_{2},\Sigma_{2})}\right)$
$\displaystyle=\frac{1}{2}\big{[}\operatorname{tr}(\Sigma_{2}^{-1}\Sigma_{1}-I)+(\mu_{2}-\mu_{1})^{\top}\Sigma_{2}^{-1}(\mu_{2}-\mu_{1})-\ln\det(\Sigma_{2}^{-1}\Sigma_{1})\big{]}.$

Moreover, suppose that all the eigenvalues of $\Sigma_{2}^{-1}\Sigma_{1}$ are
at least $\frac{1}{2}$. Then

$\displaystyle
D_{\mathrm{KL}}\left({\mathcal{N}(\mu_{1},\Sigma_{1})}\parallel{\mathcal{N}(\mu_{2},\Sigma_{2})}\right)$
$\displaystyle\leq\frac{1}{2}\big{[}\|\Sigma_{2}^{-1/2}\Sigma_{1}\Sigma_{2}^{-1/2}-I\|_{F}^{2}+(\mu_{2}-\mu_{1})^{\top}\Sigma_{2}^{-1}(\mu_{2}-\mu_{1})\big{]}$

The first assertion of 2.1 is a well-known formula for the KL divergence
between Normal distributions. The proof of the second assertion can be found
in Appendix B.1.

######  Lemma 2.2 (Pinsker’s Inequality [42, Lemma 2.5]).

For any two distributions $\mathcal{D}_{1},\mathcal{D}_{2}$,
$d_{\mathrm{TV}}\left({\mathcal{D}_{1}},{\mathcal{D}_{2}}\right)^{2}\leq\frac{1}{2}D_{\mathrm{KL}}\left({\mathcal{D}_{1}}\parallel{\mathcal{D}_{2}}\right)$.

The following lemma is a standard bound on the (upper) tails of a Gaussian.

######  Lemma 2.3 (e.g. [44, Proposition 2.1.2]).

Let $g\sim\mathcal{N}(0,\sigma^{2})$. For all $x>0$,
$\mathbb{P}\left[g\geq\sigma x\right]\leq e^{-x^{2}/2}$.

We require the following concentration inequality which can be easily derived
from well-known techniques (see [45, §2.1.3]). For completeness, we provide a
proof in Appendix B.2.

######  Lemma 2.4.

Let $g\sim\mathcal{N}(0,I_{d})$ and $A$ be a symmetric matrix. Then for all
$x>0$,

$\mathbb{P}\left[g^{\top}Ag\geq\operatorname{tr}(A)+4\|A\|_{F}\cdot\sqrt{x}+4\|A\|\cdot
x\right]\leq\exp(-x).$

###  2.2 Differential Privacy

In this section, we let $\mathcal{X},\mathcal{Y}$ denote sets and
$n\in\mathbb{N}$. Two datasets
$D=(X_{1},\ldots,X_{n}),D^{\prime}=(X_{1},\ldots,X_{n})\in\mathcal{X}^{n}$ are
said to be _neighbouring_ if $d_{H}(D,D^{\prime})\leq 1$ where $d_{H}$ denotes
Hamming distance, i.e. $d_{H}(D,D^{\prime})=|\\{i\in[d]\,:\,X_{i}\neq
X_{i}^{\prime}\\}|$.

######  Definition 2.5 ([19, 17]).

A mechanism $\mathcal{M}\colon\mathcal{X}^{n}\to\mathcal{Y}$ is said to be
$(\varepsilon,\delta)$-differentially private if for all neighbouring datasets
$D,D^{\prime}\in\mathcal{X}^{n}$ and all $Y\subseteq\mathcal{Y}$, we have

$\mathbb{P}\left[\mathcal{M}(D)\in Y\right]\leq
e^{\varepsilon}\cdot\mathbb{P}\left[\mathcal{M}(D^{\prime})\in
Y\right]+\delta.$

In this paper, we use DP as shorthand for differentially private or
differential privacy, depending on context.

######  Remark 2.6.

The astute reading may notice that the algorithms we present in this paper are
_not_ invariant under permutation of the input. Given this, it may seem more
desirable to define differential privacy where we view datasets as unordered
sets instead of ordered tuples. However, one can show that if one has an
$(\varepsilon,\delta)$-DP algorithm for ordered tuples than one also has a
$(\varepsilon,\delta)$-DP for unordered sets by randomly permuting the input
[6, Proposition D.6].

######  Definition 2.7.

Let $\mathcal{D}_{1},\mathcal{D}_{2}$ be two continuous distributions defined
on $\mathbb{R}^{d}$ and let $f_{1},f_{2}$ be the respective density functions.
We use
$\mathcal{L}_{{\mathcal{D}_{1}}\parallel{\mathcal{D}_{2}}}\colon\mathbb{R}^{d}\to\mathbb{R}$
to denote the logarithm of the likelihood ratio, i.e. for any
$x\in\mathbb{R}^{d}$,

$\mathcal{L}_{{\mathcal{D}_{1}}\parallel{\mathcal{D}_{2}}}(x)\coloneqq\ln\frac{f_{1}(x)}{f_{2}(x)}.$
(2.1)

######  Remark 2.8.

Note that
$\mathbb{E}_{Y\sim\mathcal{D}_{1}}\left[\mathcal{L}_{{\mathcal{D}_{1}}\parallel{\mathcal{D}_{2}}}(Y)\right]=D_{\mathrm{KL}}\left({\mathcal{D}_{1}}\parallel{\mathcal{D}_{2}}\right)$.

######  Definition 2.9.

Let $\mathcal{D}_{1},\mathcal{D}_{2}$ be two distributions defined on a set
$\mathcal{X}$. Then $\mathcal{D}_{1},\mathcal{D}_{2}$ are said to be
$(\varepsilon,\delta)$-indistinguishable if for all measurable
$S\subseteq\mathcal{X}$, we have

$\mathbb{P}_{Y\sim\mathcal{D}_{1}}\left[Y\in S\right]\leq
e^{\varepsilon}\mathbb{P}_{Y\sim\mathcal{D}_{2}}\left[Y\in
S\right]+\delta\quad\text{and}\quad\mathbb{P}_{Y\sim\mathcal{D}_{2}}\left[Y\in
S\right]\leq e^{\varepsilon}\mathbb{P}_{Y\sim\mathcal{D}_{1}}\left[Y\in
S\right]+\delta.$

In light of Definition 2.9, a mechanism $\mathcal{M}$ is
$(\varepsilon,\delta)$-DP if and only if for every neighbouring datasets
$D,D^{\prime}$, the distributions given by $\mathcal{M}(D)$ and
$\mathcal{M}(D^{\prime})$ are $(\varepsilon,\delta)$-indistinguishable.

In this paper, we make use of the following sufficient condition for
indistinguishability; a proof can be found in Appendix B.3.

######  Lemma 2.10.

Let $\mathcal{D}_{1},\mathcal{D}_{2}$ be continuous distributions defined on
$\mathbb{R}^{d}$. If

$\mathbb{P}_{Y\sim\mathcal{D}_{1}}\left[\mathcal{L}_{{\mathcal{D}_{1}}\parallel{\mathcal{D}_{2}}}(Y)\geq\varepsilon\right]\leq\delta\quad\text{and}\quad\mathbb{P}_{Y\sim\mathcal{D}_{2}}\left[\mathcal{L}_{{\mathcal{D}_{2}}\parallel{\mathcal{D}_{1}}}(Y)\geq\varepsilon\right]\leq\delta$

then $\mathcal{D}_{1},\mathcal{D}_{2}$ are
$(\varepsilon,\delta)$-indistinguishable.

The following claim follows from the probability density function of Normal
distributions.

######  Claim 2.11.

Let $\mu_{1},\mu_{2}\in\mathbb{R}^{d}$ and $\Sigma_{1},\Sigma_{2}\succ 0$.
Then

$\mathcal{L}_{{\mathcal{N}(\mu_{1},\Sigma_{1})}\parallel{\mathcal{N}(\mu_{2},\Sigma_{2})}}(x)=\ln\left(\frac{\det(\Sigma_{1})^{-1/2}\exp\left(-\frac{1}{2}(x-\mu_{1})^{\top}\Sigma_{1}^{-1}(x-\mu_{1})\right)}{\det(\Sigma_{2})^{-1/2}\exp\left(-\frac{1}{2}(x-\mu_{2})^{\top}\Sigma_{2}^{-1}(x-\mu_{2})\right)}\right).$
(2.2)

In this paper, we make use of the truncated Laplace distribution, denoted by
$\operatorname{TLap}(\Delta,\varepsilon,\delta)$ whose probability density
function is given by

$f_{\operatorname{TLap}(\Delta,\varepsilon,\delta)}(x)\coloneqq\begin{cases}Be^{-|x|/\lambda}&x\in[-A,A]\\\
0&x\notin[-A,A]\end{cases},$

where $\lambda=\frac{\Delta}{\varepsilon}$,
$A=\frac{\Delta}{\varepsilon}\ln\left(1+\frac{e^{\varepsilon}-1}{2\delta}\right)$,
$B=\frac{1}{2\lambda(1-e^{-A/\lambda})}$.

######  Theorem 2.12 ([22, Theorem 1]).

Suppose that $q\colon\mathcal{X}\to\mathbb{R}$ is a function with sensitivity
$\Delta$. Then the mechanism $q(x)+Y$ where
$Y\sim\operatorname{TLap}(\Delta,\varepsilon,\delta)$ is
$(\varepsilon,\delta)$-DP.

##  3 A framework for private statistical estimation

In this section, we describe an efficient and fairly general reduction from
private statistical estimation to the non-private counterpart. Our reduction
is inspired by the Propose-Test-Release framework [18] and the Subsample-And-
Aggregate framework [36]. At a very high level, we partition the data set into
multiple disjoint subsets and run the non-private estimation method on each
subset. We then privately check whether most of these solutions are “close” to
each-other. If this is not the case then the method fails. Otherwise, we
compute a weighted average of the “concentrated” solutions. Roughly speaking,
solutions which are “close” to many other candidate solutions (say $0.7$
fraction of them) are given weight $1$ while solutions that are “close” to
fewer candidate solutions (say less than $0.6$ fraction of them) are given
weight $0$. To help with sensitivity, we then linearly interpolate the weight
for solutions that are in the middle. Finally, we output a “noisy/masked”
version of this average. In order for the reduction to work in rather general
scenarios (such as learning the high dimensional Gaussians with respect the TV
distance), we first define the notion of convex semimetric spaces. This helps
to formalize the notion of “closeness” between solutions as well as the
weighted average over a number of solutions.

###  3.1 Convex semimetric spaces and masking mechanisms

In this paper, we work with the following notion of a convex semimetric space.
The key property to keep in mind is that for semimetric spaces, we only have
an _approximate_ triangle inequality, as long as the points are significantly
close together.

######  Definition 3.1 (Convex semimetric space).

Let $\mathcal{Y}$ be a convex set and let
$\operatorname{dist}\colon\mathcal{Y}\times\mathcal{Y}\to\mathbb{R}_{\geq 0}$.
We say $(\mathcal{Y},\operatorname{dist})$ is a convex semimetric space if
there exists absolute constants $t\geq 1$, $\phi\geq 0$, and $r>0$ such that
for every $k\in\mathbb{N}$ and every
$Y,Y_{1},Y_{2},\ldots,Y_{k}\in\mathcal{Y}$, the following conditions hold.

  1. 1.

$\operatorname{dist}(Y,Y)=0$ and $\operatorname{dist}(Y_{1},Y_{2})\geq 0$.

  2. 2.

Symmetry. $\operatorname{dist}(Y_{1},Y_{2})=\operatorname{dist}(Y_{2},Y_{1})$.

  3. 3.

$t$-approximate $r$-restricted triangle inequality. If
$\operatorname{dist}(Y_{1},Y_{2}),\operatorname{dist}(Y_{2},Y_{3})\leq r$ then
$\operatorname{dist}(Y_{1},Y_{3})\leq
t\cdot(\operatorname{dist}(Y_{1},Y_{2})+\operatorname{dist}(Y_{2},Y_{3}))$.

  4. 4.

Convexity. For all $\alpha\in\Delta_{k}$,
$\operatorname{dist}\left(\sum\alpha_{i}Y_{i},Y\right)\leq\sum\alpha_{i}\operatorname{dist}(Y_{i},Y)$.

  5. 5.

$\phi$-Locality. For all $\alpha,\alpha^{\prime}\in\Delta_{k}$,
${\operatorname{dist}\left(\sum_{i}{\alpha_{i}}Y_{i},\sum_{i}{\alpha_{i}^{\prime}Y_{i}}\right)\leq\sum_{i}|\alpha_{i}-\alpha_{i}^{\prime}|(\phi+\max_{i,j}\left(\operatorname{dist}(Y_{i},Y_{j})\right)}$

where $\Delta_{k}$ denotes the $k$-dimensional probability simplex. When $r$
is unspecified we take it to mean that $r=\infty$ and refer to it as a
$t$-approximate triangle inequality.

The following technical lemma (whose proof appears in Appendix C) is helpful
for learning covariance matrices.

######  Lemma 3.2.

Let $\mathcal{S}^{d}$ be the set of all $d\times d$ positive definite
matrices. For $A,B\in\mathcal{S}^{d}$ let
$\operatorname{dist}(A,B)=\max\\{\|A^{-1/2}BA^{-1/2}-I\|,\|B^{-1/2}AB^{-1/2}-I\|\\}$.
Then $(\mathcal{S}^{d},\operatorname{dist})$ is a convex semimetric which
satisfies a $(3/2)$-approximate $1$-restricted triangle inequality and
$1$-locality.

######  Definition 3.3 (Masking mechanism).

Let $(\mathcal{Y},\operatorname{dist})$ be a convex semimetric space. A
randomized function $\mathcal{B}\colon\mathcal{Y}\to\mathcal{Y}$ is a
$(\gamma,\varepsilon,\delta)$-masking mechanism for
$(\mathcal{Y},\operatorname{dist})$ if for all $Y_{1},Y_{2}\in\mathcal{Y}$
satisfying $\operatorname{dist}(Y_{1},Y_{2})\leq\gamma$, we have that
$\mathcal{B}(Y_{1}),\mathcal{B}(Y_{2})$ are
$(\varepsilon,\delta)$-indistinguishable.

######  Remark 3.4.

It is straightforward to show that for $p\geq 1$, the vector space
$\mathbb{R}^{d}$ equipped with
$\operatorname{dist}(Y_{1},Y_{2})=\|Y_{1}-Y_{2}\|_{p}$ forms a convex
semimetric space (see Proposition 7.1). In particular, the Laplace mechanism
[19] is a masking mechanism for $(\mathbb{R}^{d},\ell_{1})$ and the Gaussian
mechanism is a masking mechanism for $(\mathbb{R}^{d},\ell_{2})$. In Section
5, we design a masking mechanism for the set of covariance matrices with
respect to
$\operatorname{dist}(A,B)=\max\\{\|A^{-1/2}BA^{-1/2}-I\|,\|B^{-1/2}AB^{-1/2}-I\|\\}$.

######  Definition 3.5 (Concentration of a masking mechanism).

Let $(\mathcal{Y},\operatorname{dist})$ be a convex semimetric space and
$\mathcal{B}$ be $(\gamma,\varepsilon,\delta)$-masking mechanism for it. We
say $\mathcal{B}$ is $(\alpha,\beta)$ concentrated if for all
$Y\in\mathcal{Y}$,
$\mathbb{P}\left[\operatorname{dist}(\mathcal{B}(Y),Y)>\alpha\right]\leq\beta$.

###  3.2 The Private Populous Mean Estimator

In this section we outline an efficient and fairly general reduction from
private statistical estimation to the non-private counterpart. We use
$\mathcal{A}\colon\mathcal{X}^{*}\to\mathcal{Y}$ to denote a generic
(efficient) non-private mechanism which takes a collection of samples as
inputs and outputs an element in $\mathcal{Y}$. We assume that
$(\mathcal{Y},\operatorname{dist})$ is a convex semimetric space, and we have
access to (an efficient) masking mechanism for
$(\mathcal{Y},\operatorname{dist})$. Intuitively, the error of the estimator
is measured by $\operatorname{dist}(\hat{Y},Y^{*})$, where $\hat{Y}$ is our
estimation while $Y^{*}$ is the ground truth.

At a high-level, the algorithm works as follows. We split the dataset into $k$
disjoint groups and run the non-private algorithm $\mathcal{A}$ on each group
$i$ to obtain an output $Y_{i}$. We then consider a ball of radius $r/t$
around $Y_{i}$ with respect to $\operatorname{dist}$444$r$ and $t$ are chosen
in a way that $(\mathcal{Y},\operatorname{dist})$ satisfies $r$-restricted
$t$-approximate triangle inequality. and compute $q_{i}$, the fraction of
$Y_{j}$’s that are inside this ball. We then look at the average value of the
$q_{i}$, i.e., $Q=\frac{1}{k}\sum_{i\in[k]}q_{i}$. Intuitively, $Q$
efficiently measures the “stability” of $\mathcal{A}$ across the data set
partitions. Namely, $Q$ is large when most of the $Y_{i}$’s are concentrated
in a small “popular region”, and small when $Y_{i}$’s are scattered around
$\mathcal{Y}$. If $Q$ is small, then $\mathcal{A}$ is not sufficiently stable
and we can announce failure privately using a noisy threshold. If $Q$ is
large, then we take a weighted average of $Y_{i}$’s that are in the popular
region and then output a “masked” version of this average. The following lemma
establishes the privacy and accuracy of Algorithm 1 and its proof can be found
in Appendix C.2.

Algorithm 1 Private to non-private reduction via populous mean estimation

Input: Dataset $D=(X_{1},\ldots,X_{m})$; parameters
$\varepsilon,\delta,r,\phi>0$; $k\leq m$; $t\geq 2$.

1:Let $s\leftarrow\lfloor m/k\rfloor$.

2:For $i\in[k]$, let
$Y_{i}\leftarrow\mathcal{A}(\\{X_{\ell}\\}_{\ell=(i-1)s+1}^{is})$.

3:For $i\in[k]$, let
$q_{i}\leftarrow\frac{1}{k}|\\{j\in[k]\,:\,\operatorname{dist}(Y_{i},Y_{j})\leq
r/t\\}|$.

4:Let $Q\leftarrow\frac{1}{k}\sum_{i\in[k]}q_{i}$.

5:Let $Z\sim\operatorname{TLap}(2/k,\varepsilon,\delta)$.

6:Let $\widehat{Q}\leftarrow Q+Z$.

7:If
$\widehat{Q}<0.8+\frac{2}{k\varepsilon}\ln\left(1+\frac{e^{\varepsilon}-1}{2\delta}\right)$,
fail and return $\perp$.

8:For $i\in[k]$, let $w_{i}=\min\left(1,10\max(0,q_{i}-0.6)\right)$

9:$\mu=(\sum_{i=1}^{k}w_{i}Y_{i})/(\sum_{i=1}^{k}w_{i})$

10:Return $\mathcal{B}(\mu)$.

######  Lemma 3.6.

Let $(\mathcal{Y},\operatorname{dist})$ be a convex semimetric space that
satisfies $r$-restricted $t$-approximate triangle inequality and
$\phi$-locality for some $r>0,t\geq 1,\phi\geq 0$. Let $\mathcal{B}$ be a
$(400(r+\phi)/k,\varepsilon,\delta)$-masking mechanism with respect to
$(\mathcal{Y},\operatorname{dist})$. Then

  1. 1.

Privacy. For $k\geq 140$, Algorithm 1 is
$(2\varepsilon,4e^{\varepsilon}\delta)$-DP.

  2. 2.

Utility. Assume
$k\geq\frac{20}{\varepsilon}\ln\left(1+\frac{e^{\varepsilon}-1}{2\delta}\right)$.
For $\alpha_{1}\leq r/t$ assume that there exist $Y^{*}\in\mathcal{Y}$ such
that $\operatorname{dist}(Y_{i},Y^{*})\leq\alpha_{1}/t$ for all $i\in[k]$. If
$\mathcal{B}$ is ($\alpha_{2}/t,\beta$)-concentrated then the output $\hat{Y}$
of Algorithm 1 satisfies

$\mathbb{P}\left[\operatorname{dist}(\hat{Y},Y^{*})>\alpha_{1}+\alpha_{2}\right]\leq\beta$

######  Remark 3.7.

For privacy, 3.6 requires that $k$ is larger than some constant. If this is
not the case, we can always add “dummy” points into the dataset.

######  Remark 3.8.

The computational complexity of Algorithm 1 is dominated by these
computations:

  * •

$k$ queries to the oracle $\mathcal{A}$ (with inputs of size $m/k$).

  * •

$k^{2}$ computations of the $\operatorname{dist}(,.,)$ function.

  * •

One query to the masking mechanism $\mathcal{B}$.

In particular, if $\mathcal{A}$, $\mathcal{B}$, and the algorithm for
computing $\operatorname{dist}$ all run in polynomial time then Algorithm 1
runs in polynomial time.

##  4 Learning the subspace

As a simple application of our technique, we consider the problem of privately
learning a subspace. More specifically, if the samples are drawn from
$\mathcal{N}(0,\Sigma)$ then the goal is to learn the range of $\Sigma$.
Recently, Singhal and Steinke [39] also studied the problem of privately
learning a subspace albeit under slightly different assumption than ours. For
learning a subspace _exactly_ , they gave an algorithm that was sample optimal
but not computationally efficient. For learning a subspace _approximately_ ,
they provided an algorithm that was sample efficient and computationally
efficient. The latter is used (with some modifications) in [25] to efficiently
learn the covariance matrix of a Gaussian distribution. In this section, we
give a computationally-efficient algorithm for learning a subspace _exactly_
with a slightly different sample complexity (albeit with different
assumptions). Formally, we prove the following theorem.

######  Theorem 4.1.

There exists a mechanism $\mathcal{M}$ satisfying the following. For all
$m,d\in\mathbb{N}$, $\varepsilon,\delta\in(0,1)$, and dataset
$D=(X_{1},\ldots,X_{m})$:

  1. 1.

Given $\varepsilon,\delta,D$ as input, $\mathcal{M}$ is
$(2\varepsilon,4e^{\varepsilon}\delta)$-DP and runs in time $O(md)$.

  2. 2.

If $\Sigma\succeq 0$,
$X_{1},\ldots,X_{m}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}\mathcal{N}(0,\Sigma)$,
and
$m\geq\max\\{140,\lceil\frac{20}{\varepsilon}\ln(1+\frac{e^{\varepsilon}-1}{2\delta})\rceil\\}\cdot
d=O(d\log(1/\delta)/\varepsilon)$ then, almost surely, $\mathcal{M}$ outputs
the orthogonal projection matrix onto the range of $\Sigma$.

The mechanism that achieves Theorem 4.1 appears in Algorithm 2 and is designed
using the template given in Algorithm 1. In this section, we define

$\operatorname{dist}(A,B)=\begin{cases}0&\text{if $A=B$ are orthogonal
projection matrices}\\\ \infty&\text{otherwise}\end{cases},$

where $A,B$ are two matrices in $\mathbb{R}^{d\times d}$. Note that this makes
the choice of $r,t,\phi$ moot since the distance between any two matrices is
either $0$ or $\infty$. Finally, recall that orthogonal projection matrices
are unique. In other words, if $P_{1}$ and $P_{2}$ are both orthogonal
projection matrices onto a subspace $S$ then $P_{1}=P_{2}$.

Algorithm 2 A private algorithm for learning the subspace.

Input: Dataset $D=(X_{1},\ldots,X_{m})$; parameters
$\varepsilon,\delta\in(0,1)$.

1:$\triangleright$ Parameter settings.

2:$r\leftarrow 1$, $t\leftarrow 1$, $\phi\leftarrow 1$. $\triangleright$ Not
used in Algorithm 1.

3:$k\leftarrow\max\\{30,\frac{20}{\varepsilon}\ln\left(1+\frac{e^{\varepsilon}-1}{2\delta}\right)\\}$

4:function $\mathcal{A}$($X_{1},\ldots,X_{s}$)

5: Let $P$ be orthogonal projection matrix onto
$\operatorname{span}\\{X_{1},\ldots,X_{s}\\}$.

6: Return $P$.

7:end function

8:function $\mathcal{B}$($P$)

9: Return $P$.

10:end function

###### Proof of Theorem 4.1.

Using the convention $\infty\leq\infty$, it is easy to check that the
$\operatorname{dist}$ function we defined above forms a convex semimetric
space for matrices (see Definition 3.1). Also, the identity function is a
$(\gamma,\varepsilon,\delta)$-masking mechanism for all
$\gamma,\varepsilon,\delta\in(0,1)$ since for all $A,B$,
$\operatorname{dist}(A,B)<\infty$ only if $A,B$ are both orthogonal projection
matrices and $A=B$. Therefore, the privacy statement follows from 3.6.

For the utility, note that the masking mechanism is the identity function, so
it is $(0,0)$-concentrated. Since we set
$k=\max\left\\{140,\left\lceil\frac{20}{\varepsilon}\ln\left(1+\frac{e^{\varepsilon}-1}{2\delta}\right)\right\rceil\right\\}$,
our assumption on $m$ implies that $m\geq kd$. In particular, $s\geq d$ (see
Algorithm 1). Therefore, with probability $1$, the output of
$\mathcal{A}(\\{X_{\ell}\\}_{\ell=(i-1)s+1}^{is})$ is the orthogonal
projection matrix onto the range of $\Sigma$ for all $i\in[k]$ almost surely.
Therefore, the utility conditions of 3.6 are satisfied with
$\alpha_{1}=\alpha_{2}=\beta=0$. ∎

##  5 An efficient algorithm for learning the covariance matrix

In this section we show how to use the general reduction to learn high-
dimensional Gaussian distributions with respect to the total variation
distance. Here, the main hurdle is to find an approximation to the covariance
matrix for the zero-mean case. Based on 3.2, the following
$\operatorname{dist}$ function forms a semimetric space for positive definite
matrices.

$\operatorname{dist}(\Sigma_{1},\Sigma_{2})=\begin{cases}\max\left(\|\Sigma_{2}^{-1/2}\Sigma_{1}\Sigma_{2}^{-1/2}-I_{d}\|,\|\Sigma_{1}^{-1/2}\Sigma_{2}\Sigma_{1}^{-1/2}-I_{d}\|\right)&\text{if
$\operatorname{rank}\Sigma_{1}=\operatorname{rank}\Sigma_{2}=d$}\\\
\infty&\text{otherwise}\end{cases}.$ (5.1)

The main steps to use the reduction are to define a masking mechanism and show
that it is concentrated. We start by stating the main algorithmic idea of this
section, which is a masking mechanism with respect to the
$\operatorname{dist}$ function in Eq. (5.1).

Before we state our masking mechanism, we provide some intuition on its
design. First, imagine that we were working on $\mathbb{R}^{d^{2}}$ with the
standard Euclidean metric. Let $G\sim\mathcal{N}(0,I_{d^{2}})$ and consider
the Gaussian mechanism $\mathcal{B}(X)=X+\eta G$. If
$X_{1},X_{2}\in\mathbb{R}^{d^{2}}$ and
$\|X_{1}-X_{2}\|_{2}^{2}\leq\gamma^{2}d$, where
$\gamma\lesssim\varepsilon\eta/\sqrt{d\ln(1/\delta)}$, then a standard
analysis shows that $\mathcal{B}(X_{1}),\mathcal{B}(X_{2})$ are
$(\varepsilon,\delta)$-indistinguishable.

The masking mechanism in 5.1 is attempting to mimic the Gaussian mechanism
where we view $d\times d$ matrices as vectors of length $d^{2}$. A first
attempt may be to take a matrix $\Sigma$, compute
$\widehat{\Sigma}=\Sigma^{1/2}+\eta G$, where $G$ is a $d\times d$ matrix with
$\mathcal{N}(0,1)$ entries, and return
$\widehat{\Sigma}\widehat{\Sigma}^{\top}$. However, this does not quite work
for covariance estimation since we require the noise to have a difference
scale in different directions, i.e. the noise along each eigenvector should
scale with the corresponding eigenvalue. Thus, we use an empirically scaled
Gaussian mechanism where the noise we add is shaped and scaled by the input
itself.

######  Lemma 5.1.

For a matrix $\Sigma\succ 0$, define $\mathcal{B}(\Sigma)=\Sigma^{1/2}(I+\eta
G)(I+\eta G)^{\top}\Sigma^{1/2}$ where $G\in\mathbb{R}^{d\times d}$ is a
$d\times d$ matrix with independent $\mathcal{N}(0,1)$ entries.

Then for every $\eta>0$ and $\varepsilon,\delta\in(0,1]$, $\mathcal{B}$ is a
($\gamma,\varepsilon,\delta$)-masking mechanism with respect to
$\operatorname{dist}$ for

$\gamma=\min\left\\{\sqrt{\frac{\varepsilon}{2d(d+1/\eta^{2})}},\frac{\varepsilon}{8d\sqrt{\ln(2/\delta)}},\frac{\varepsilon}{8\ln(2/\delta)},\frac{\varepsilon\eta}{12\sqrt{d}\sqrt{\ln(2/\delta)}}\right\\}.$

The proof of 5.1 can be found in Appendix D.1. Next, we show that the masking
mechanism in 5.1 is concentrated.

######  Lemma 5.2.

Let $\Sigma\succ 0$ and set
$\eta=\frac{1}{C_{1}(\sqrt{d}+\sqrt{\ln(4/\beta)})}$ for a sufficiently large
(universal) constant $C_{1}>0$. Then the masking mechanism
$\mathcal{B}(\Sigma)$ defined in 5.1 is $(1/100,\beta/2)$-concentrated.

###### Proof.

Let $\widehat{\Sigma}=\mathcal{B}(\Sigma)=\Sigma^{1/2}(I+\eta G)(I+\eta
G)^{\top}\Sigma^{1/2}$. First, we prove that

$\|\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2}-I_{d}\|\leq 1/400,$

with probability $1-\beta/2$. Indeed, we have that, for a sufficiently large
constant $C^{\prime}>0$,

$\displaystyle\|\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2}-I_{d}\|$
$\displaystyle=\|(I_{d}+\eta G)(I+\eta G)^{\top}-I_{d}\|$ $\displaystyle\leq
2\eta G+\eta^{2}\|GG^{\top}\|$ $\displaystyle\leq 2\eta
C^{\prime}(\sqrt{d}+\sqrt{\ln(4/\beta)})+\eta^{2}(C^{\prime}(\sqrt{d}+\sqrt{\ln(4/\beta)}))^{2},$

where we used the fact that $\|G\|\leq C(\sqrt{d}+\sqrt{\ln(4/\beta)})$ with
probability $1-\beta/2$ (see Theorem A.5). Taking
$\eta\leq\frac{1}{1600C(\sqrt{d}+\sqrt{\ln(4/\beta)})}$ gives that

$\|\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2}-I_{d}\|\leq 1/400,$

as desired.

In order to show that $\operatorname{dist}(\Sigma,\widehat{\Sigma})$, it
remains to show that

$\|\widehat{\Sigma}^{-1/2}\Sigma\widehat{\Sigma}^{-1/2}-I_{d}\|\leq 1/100.$

This follows directly from A.6. ∎

Algorithm 3 summarizes the modules needed for 3.6 to work, including the
definitions of the mechanism $\mathcal{A},\mathcal{B}$ as well as the required
parameters. Note that $\mathcal{A}$ returns just the empirical covariance
matrix (although we could easily replace it by any algorithm which
approximates the covariance matrix). One last piece is to show that the output
of the non-private $\mathcal{A}$ is accurate enough.

######  Lemma 5.3.

There is a universal constant $C_{2}$ such that the following holds. Suppose
$s\geq C_{2}(d+\ln(4/\beta))$. Let
$X_{1},\ldots,X_{s}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}\mathcal{N}(0,\Sigma)$
and $\widehat{\Sigma}=\sum_{\ell=1}^{s}X_{\ell}X_{\ell}^{\top}$. Then
$\operatorname{dist}(\Sigma,\widehat{\Sigma})\leq 1/100$ with probability
$1-\beta/2$.

###### Proof.

Applying a standard concentration inequality (see A.4), we have
$\|\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2}-I_{d}\|\leq
C\left(\sqrt{\frac{d+\ln(4/\beta)}{s}}+\frac{d+\ln(4/\beta)}{s}\right)$.
Choosing $s\geq 800^{2}C^{2}(d+\ln(4/\beta))$ gives that
$\|\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2}-I_{d}\|\leq 1/400$. Applying A.6
gives that $\operatorname{dist}(\widehat{\Sigma},\Sigma)\leq 1/100$ with
probability $1-\beta/2$ as required. ∎

Algorithm 3 A private algorithm for learning the covariance matrix up to a
constant factor with respect to the spectral distance

Input: Dataset $D=(X_{1},\ldots,X_{m})$; parameters
$\varepsilon,\delta,\alpha,\beta\in(0,1)$

1:$\triangleright$ Some parameter settings.

2:$\eta\leftarrow\frac{1}{C_{1}(\sqrt{d}+\sqrt{\ln(4/\beta)})}$
$\triangleright$ $C_{1}$ from 5.2

3:$\gamma\leftarrow\min\left\\{\sqrt{\frac{\varepsilon}{2d(d+1/\eta^{2})}},\frac{\varepsilon}{8d\sqrt{\ln(2/\delta)}},\frac{\varepsilon}{8\ln(2/\delta)},\frac{\varepsilon\eta}{12\sqrt{d\ln(2/\delta)}}\right\\}$.

4:$r\leftarrow 1$, $\phi\leftarrow 1$

5:$k\leftarrow\max\left\\{\frac{400(r+\phi)}{\gamma},\frac{20}{\varepsilon}\ln\left(1+\frac{e^{\varepsilon}-1}{2\delta}\right)\right\\}$

6:$t\leftarrow 3/2$

7:function $\mathcal{A}$($X_{1},\ldots,X_{s}$)

8: $\widehat{\Sigma}=\frac{1}{s}\sum_{\ell=1}^{s}X_{\ell}X_{\ell}^{\top}$

9: Return $\widehat{\Sigma}$.

10:end function

11:function $\mathcal{B}$($\widehat{\Sigma}$)

12: Let $G$ be $d\times d$ matrix with independent $\mathcal{N}(0,1)$ entries.

13: Return $\widehat{\Sigma}^{1/2}(I+\eta G)(I+\eta
G)^{\top}\widehat{\Sigma}^{1/2}$.

14:end function

Now we are ready to apply 3.6. The following theorem shows that we can learn
the covariance matrix up to a constant distance with respect to the spectral
distance using more or less $O(d^{2}/\varepsilon)$ samples.

######  Theorem 5.4.

Applying 3.6 with parameters specified in Algorithm 3 gives a mechanism
$\mathcal{M}$ satisfying the following. For all $m,d\in\mathbb{N}$,
$\varepsilon,\delta,\beta\in(0,1)$, and dataset $D=(X_{1},\ldots,X_{m})$:

  1. 1.

Given $\varepsilon,\delta,\beta,D$ as input, $\mathcal{M}$ is
$(2\varepsilon,4e^{\varepsilon}\delta)$-DP and runs in time
$O(md^{2}+kd^{3})$.

  2. 2.

If $\Sigma\succ 0$,
$X_{1},\ldots,X_{m}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}\mathcal{N}(0,\Sigma)$,
and
$m=\widetilde{\Omega}\left(\frac{d^{2}\ln(1/\beta)^{3/2}\ln(1/\delta)^{1/2}}{\varepsilon}\right)$
then with probability at least $1-\beta$, $\mathcal{M}$ outputs a matrix
$\widehat{\Sigma}$ such that $\operatorname{dist}(\Sigma,\widehat{\Sigma})\leq
1/10$.

###### Proof.

3.2 shows that $(\mathcal{S}^{d},\operatorname{dist})$ satisfies the
assumptions of Definition 3.1. 5.1 (with $r=\phi=1,t=3/2$) shows that
$\mathcal{B}$ (as defined in 5.1 and Algorithm 3) as a
$(\gamma,\varepsilon,\delta)$-masking mechanism. Since $k\geq
400(r+\phi)/\gamma$ (so $400(r+\phi)/k\leq\gamma$) it follows that
$\mathcal{B}$ is also a $(400(r+\phi)/k,\varepsilon,\delta)$-masking
mechanism. Therefore, privacy immediately follows from the first assertion of
3.2.

Next, if $m\geq C_{2}(d+\ln(4k/\beta))\cdot
k=\widetilde{O}\left(\frac{d(d\sqrt{\ln(1/\delta)}+\ln(1/\delta))\ln(1/\beta)^{3/2}}{\varepsilon}\right)$
then $\widehat{\Sigma}_{i}=\mathcal{A}(\\{X_{\ell}\\}^{is}_{\ell=(i-1)s+1})$
satisfies $\operatorname{dist}(\Sigma,\widehat{\Sigma}_{i})\leq 1/10\leq r/t$
for all $i\in[k]$ with probability $1-\beta/2$ by 5.3 and a union bound over
$k$. In addition, $\mathcal{B}$ is $(1/100,\beta/2)$-concentrated.

To conclude, let $Y_{i}=\widehat{\Sigma}_{i}$, $Y^{*}=\Sigma$, and
$\widehat{Y}=\widehat{\Sigma}$ be the output of Algorithm 1 with the modules
specified in Algorithm 3. Then, taking $\alpha_{1}=\alpha_{2}=1/25$ in 3.2
gives that $\operatorname{dist}(\widehat{\Sigma},\Sigma)\leq 1/10$ with
probability $1-\beta$. ∎

Assuming $\Sigma\succ 0$, Theorem 5.4 gives us a covariance matrix
$\widehat{\Sigma}$ such that
$0.9I_{d}\preceq\widehat{\Sigma}^{-1/2}\Sigma\widehat{\Sigma}^{-1/2}\preceq
1.1I_{d}$. Note that if $X\sim\mathcal{N}(0,\Sigma)$ then
$\widehat{\Sigma}^{-1/2}X\sim\mathcal{N}(0,\widehat{\Sigma}^{-1/2}\Sigma\widehat{\Sigma}^{-1/2})$.
We can now make use of off-the-shelf private covariance methods which are
efficient when the covariance matrix is well-conditioned [5, 24]. For
concreteness, we use the following theorem which is paraphrased from [24,
Theorem 1.1], [5, Theorem 3.3].

######  Theorem 5.5 ([5, 24]).

There is a polynomial time $(\varepsilon,\delta)$-DP algorithm that given

$m=\widetilde{\Omega}\left(\left(\frac{d^{2}}{\alpha^{2}}+\frac{d^{2}\sqrt{\ln(1/\delta)}}{\alpha\varepsilon}\right)\cdot\ln(1/\beta)\right)$

samples from a Gaussian $\mathcal{N}(0,\Sigma)$ with unknown covariance matrix
$\Sigma$ such that $0.9I\preceq\Sigma\preceq 1.1I$ and outputs
$\widehat{\Sigma}$ satisfying
$\|\widehat{\Sigma}^{-1/2}\Sigma\widehat{\Sigma}^{-1/2}-I_{d}\|_{F}\leq\alpha$.
In particular,
$d_{\mathrm{TV}}\left({\mathcal{N}(0,\Sigma)},{\mathcal{N}(0,\widehat{\Sigma})}\right)\leq\alpha$.

Finally, we prove our result for general, full-rank covariance matrices. Note
that one can use Theorem 4.1 to easily extend this to the case where $\Sigma$
is not necessarily full-rank. The details can be found in Subsection 5.1.

######  Theorem 5.6.

There exists a mechanism $\mathcal{M}$ satisfying the following. For all
$m,d\in\mathbb{N}$, $\varepsilon,\delta,\alpha,\beta\in(0,1)$, and dataset
$D=(X_{1},\ldots,X_{m})$:

  1. 1.

Given $\varepsilon,\delta,\alpha,\beta,D$ as input, $\mathcal{M}$ is
$(3\varepsilon,(4e^{\varepsilon}+1)\delta)$-DP and runs in polynomial time.

  2. 2.

If $\Sigma\in\mathbb{S}^{d}$,
$X_{1},\ldots,X_{m}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}\mathcal{N}(0,\Sigma)$,
and

$m=\widetilde{\Omega}\left(\frac{(d^{2}\sqrt{\ln(1/\delta)}+d\ln(1/\delta))\ln(1/\beta)^{3/2}}{\varepsilon}+\left(\frac{d^{2}}{\alpha^{2}}+\frac{d^{2}\sqrt{\ln(1/\delta)}}{\alpha\varepsilon}\right)\cdot\ln(1/\beta)\right)$

then with probability $1-\beta$, $\mathcal{M}$ outputs a matrix
$\widehat{\Sigma}$ such that
$d_{\mathrm{TV}}\left({\mathcal{N}(0,\Sigma)},{\mathcal{N}(0,\widehat{\Sigma})}\right)\leq\alpha$.

###### Proof.

Combining Theorem 5.4 and Theorem 5.5 concludes the result. The privacy
follows from basic composition since the algorithm has only three steps. The
utility directly follows from the adopted theorems. ∎

###  5.1 An efficient algorithm for learning a Gaussian

In this subsection, we sketch the ideas needed to learn a Gaussian
$\mathcal{N}(\mu,\Sigma)$ under $(\varepsilon,\delta)$-DP with no assumptions
on $\mu,\Sigma$. Suppose we draw $2m_{1}+2m_{2}+m_{3}+m_{4}$ samples from a
Gaussian distribution $\mathcal{N}(\mu,\Sigma)$, where $m_{1}$ is the sample
complexity for learning a subspace using Algorithm 2 and $m_{2}$ is the sample
complexity for learning a full-rank covariance matrix using Algorithm 3 and
Theorem 5.5. The values of $m_{3}$ and $m_{4}$ are discussed below. We now
show how to privately and efficiently learn a Gaussian in four steps.

  1. 1.

Let $Z_{i}=\frac{X_{i}-X_{m_{1}+i}}{\sqrt{2}}$ for $i\in[m_{1}]$. We use
$Z_{1},\ldots,Z_{m_{1}}$ as input into the private subspace learning algorithm
Algorithm 2. As long as $m_{1}$ is sufficiently large, with probability
$1-\beta$, it will output a matrix $P$ which is the projection onto the column
space of $\Sigma$.

  2. 2.

For $i\in[m_{2}]$, let
$Z_{i}=\frac{X_{2m_{1}+i}-X_{2m_{1}+m_{2}+i}}{\sqrt{2}}$. Also, write
$P=U_{r}U_{r}^{\top}$ where $U_{r}$ is a $d\times r$ matrix whose rows form an
orthonormal basis for the column space of $\Sigma$ (equivalently, $P$). Note
that $U_{r}^{\top}Z_{i}\sim\mathcal{N}(0,U_{r}^{\top}\Sigma U_{r})$ and that
$U_{r}^{\top}\Sigma U_{r}$ is a full rank $r\times r$ covariance matrix. We
use Algorithm 3 to learn an $r\times r$ covariance matrix
$\widehat{\Sigma}_{r}$ such that
$\|\widehat{\Sigma}_{r}^{-1/2}U_{r}^{\top}\Sigma
U_{r}\widehat{\Sigma}_{r}^{-1/2}-I_{r}\|_{F}\leq\alpha$. In particular,
$d_{\mathrm{TV}}\left({\mathcal{N}(0,\widehat{\Sigma}_{r})},{\mathcal{N}(0,U_{r}^{\top}\Sigma
U_{r})}\right)\leq O(\alpha)$ (see Theorem 5.6).

Next, we learn the mean $\mu$ in two steps.

  1. 3.

For $i\in[m_{3}]$, let $Z_{i}=X_{2m_{1}+2m_{2}+i}$. Then
$\widehat{\Sigma}_{r}^{-1/2}U_{r}^{\top}Z_{i}\sim\mathcal{N}(\widehat{\Sigma}_{r}^{-1/2}U_{r}^{\top}\mu,\widehat{\Sigma}_{r}^{-1/2}U_{r}^{\top}\Sigma
U_{r}\widehat{\Sigma}_{r}^{-1/2})$. Note that (for $\alpha\leq 0.5$), we have
$0.5I\preceq\widehat{\Sigma}_{r}^{-1/2}U_{r}^{\top}\Sigma
U_{r}\widehat{\Sigma}_{r}^{-1/2}\preceq 1.5I$. In particular, each coordinate
of $\widehat{\Sigma}_{r}^{-1/2}U_{r}^{\top}Z_{i}$ has a variance which is
bounded below by $1/2$ and above by $1.5$. Hence, applying [29, Theorem 1.3],
we can obtain an estimate $\widehat{\mu}_{r}\in\mathbb{R}^{d}$ such that
$\|\widehat{\mu}_{r}-\widehat{\Sigma}_{r}^{-1/2}U_{r}^{\top}\mu\|_{2}\leq\alpha$
provided
$m_{3}=\widetilde{\Omega}\left(\frac{d}{\alpha^{2}}+\frac{d\ln(1/\delta)}{\alpha\varepsilon}\right)$
where $\widetilde{\Omega}$ hides
$\operatorname{polylog}(1/\beta,d,\ln(1/\delta),1/\varepsilon,1/\alpha)$. We
remark that this step can be done in polynomial time.

We now have that

$\displaystyle
d_{\mathrm{TV}}\left({\mathcal{N}(U_{r}^{\top}\mu,U_{r}^{\top}\Sigma
U_{r})},{\mathcal{N}(\widehat{\Sigma}_{r}^{1/2}\widehat{\mu}_{r},\widehat{\Sigma}_{r})}\right)$
$\displaystyle\leq
d_{\mathrm{TV}}\left({\mathcal{N}(U_{r}^{\top}\mu,U_{r}^{\top}\Sigma
U_{r})},{\mathcal{N}(U_{r}^{\top}\mu,\widehat{\Sigma}_{r})}\right)$
$\displaystyle+d_{\mathrm{TV}}\left({\mathcal{N}(U_{r}^{\top}\mu,\widehat{\Sigma}_{r})},{\mathcal{N}(\widehat{\Sigma}_{r}^{1/2}\mu_{r},\widehat{\Sigma}_{r})}\right)$
$\displaystyle=d_{\mathrm{TV}}\left({\mathcal{N}(0,U_{r}^{\top}\Sigma
U_{r})},{\mathcal{N}(0,\widehat{\Sigma}_{r})}\right)$
$\displaystyle+d_{\mathrm{TV}}\left({\mathcal{N}(\widehat{\Sigma}_{r}^{-1/2}U_{r}^{\top}\mu,I_{r})},{\mathcal{N}(\mu_{r},I_{r})}\right)$
$\displaystyle\leq O(\alpha),$

where the first line uses the triangle inequality, the second line uses the
fact that total variation distance is invariant under bijective
transformations, and the last line uses steps 2 and 3 above as well as 2.1 and
2.2 to assert that the total variation distance is $O(\alpha)$. Using the fact
that bijective transformations are invariant once again, that
$P=U_{r}U_{r}^{\top}$, and that $P\Sigma P=\Sigma$, we have

$d_{\mathrm{TV}}\left({\mathcal{N}(P\mu,\Sigma)},{\mathcal{N}(U_{r}\widehat{\Sigma}_{r}^{1/2}\widehat{\mu}_{r},U_{r}\widehat{\Sigma}_{r}U_{r}^{\top})}\right)=d_{\mathrm{TV}}\left({\mathcal{N}(U_{r}^{\top}\mu,U_{r}^{\top}\Sigma
U_{r})},{\mathcal{N}(\widehat{\Sigma}_{r}^{1/2}\widehat{\mu}_{r},\widehat{\Sigma}_{r})}\right)\leq
O(\alpha).$

All that remains is to learn $(I-P)\mu$.

  1. 4.

Now let $Z_{i}=X_{2m_{1}+2m_{2}+m_{3}+i}$ for $i\in[m_{4}]$. Note that $(I-P)$
is the orthogonal projection matrix onto the orthogonal complement of the
column space of $\Sigma$. In particular $(I-P)Z_{i}=(I-P)\mu$ with probability
$1$ if $Z_{i}\sim\mathcal{N}(0,\Sigma)$. Learning $(I-P)\mu$ can now be done
using standard techniques (for example, lazy private histograms or propose-
test-release) with sample complexity $O(\log(1/\delta)/\varepsilon)$. Let us
call the output of this step $\widehat{\mu}_{\perp}$.

With the last step (and again, using that bijective maps does not affect the
total variation distance to add in $\widehat{\mu}_{\perp}$), we have

$d_{\mathrm{TV}}\left({\mathcal{N}(\mu,\Sigma)},{\mathcal{N}(U_{r}\widehat{\Sigma}_{r}^{1/2}\widehat{\mu}_{r}+\widehat{\mu}_{\perp},U_{r}\widehat{\Sigma}_{r}U_{r}^{\top})}\right)\leq
O(\alpha)$

To conclude, we have the following theorem. We remind the reader that $\Sigma$
is not assumed to be full-rank.

######  Theorem 5.7.

There exists a mechanism $\mathcal{M}$ satisfying the following. For all
$m,d\in\mathbb{N}$, $\varepsilon,\delta,\alpha,\beta\in(0,1)$, and dataset
$D=(X_{1},\ldots,X_{m})$:

  1. 1.

Given $\varepsilon,\delta,\alpha,\beta,D$ as input, $\mathcal{M}$ is
$(O(\varepsilon),O(e^{\varepsilon}\delta))$-DP and runs in polynomial time.

  2. 2.

If $X_{1},\ldots,X_{m}\sim\mathcal{N}(\mu,\Sigma)$ and

$m=\widetilde{\Omega}\left(\frac{d^{2}}{\alpha^{2}}+\frac{d^{2}\sqrt{\ln(1/\delta)}}{\alpha\varepsilon}+\frac{d\ln(1/\delta)}{\alpha\varepsilon}\right)$

then with probability $1-\beta$, $\mathcal{M}$ outputs a matrix
$\widehat{\Sigma}\in\mathbb{R}^{d\times d}$ and a vector
$\widehat{\mu}\in\mathbb{R}^{d}$ such that
$d_{TV}(\mathcal{N}(\mu,\Sigma),\mathcal{N}(\widehat{\mu},\widehat{\Sigma}))\leq\alpha$.
Here, the $\widetilde{\Omega}$ hides
$\operatorname{polylog}(d,1/\alpha,1/\beta,1/\varepsilon,\ln(1/\delta))$.

##  6 Efficient and robust learning of the covariance matrix

From the previous sections we know that the spectral $\operatorname{dist}$
function forms a convex semimetric space for positive definite matrices.
Furthermore, we already have a masking mechanism for covariance matrices (see
5.1) which is shown to be concentrated (see 5.2). In this section we show how
one can replace the non-robust subroutines of Section 5— namely 5.3 and
Theorem 5.5—with robust ones to get a robust and efficient learning method for
Gaussians.

Following [11, 32], we use the notion of $\alpha$-corrupted sample from a
distribution $\mathcal{D}$, where an adversary receives an i.i.d. sample from
$\mathcal{D}$, alters at most an $\alpha$-fraction of them arbitrarily, and
outputs the corrupted sample. The following result is based on Theorem 4.35 in
[11] and is our main tool for robust learning of the covariance matrix.

######  Theorem 6.1 ([11]).

There are constants $C_{4}>0,0<\alpha_{0}<1/2$ such that the following holds.
There exist an efficient algorithm that receives $0<\alpha<\alpha_{0}$,
$\beta>0$, and $X_{1},\ldots,X_{m}$ as input, and outputs $\widehat{\Sigma}$
with the following guarantee: if $X_{1},\ldots,X_{m}$ are $\alpha$-corrupted
samples from $\mathcal{N}(0,\Sigma)$ for some $\Sigma\in\mathbb{S}_{>0}^{d}$,
and if
$m=\widetilde{\Omega}\left(\frac{d^{2}\log^{5}1/\beta}{\alpha^{2}}\right)$,
then $\left\|\Sigma^{1/2}\widehat{\Sigma}^{-1}\Sigma^{1/2}-I\right\|_{F}\leq
C_{4}\alpha\ln(1/\alpha)<1/9$.

Our plan is use a very similar algorithm to Algorithm 3. The key differences
are that we use Theorem 6.1 for $\mathcal{A}$ and a smaller value for $\eta$.
The proof of the following proposition is very similar to that of 5.2 and is
omitted for brevity.

######  Proposition 6.2.

Let $\Sigma\succ 0$ and $\alpha\in(0,1)$. Set
$\eta=\frac{\alpha}{C_{1}(d+\sqrt{d\ln(4/\beta)})}$ for a sufficiently large
(universal) constant $C_{1}>0$. Then the masking mechanism
$\mathcal{B}(\Sigma)$ defined in 5.1 is
$(\alpha/\sqrt{d},\beta/2)$-concentrated.

Algorithm 4 A private and robust algorithm for learning the covariance matrix

Input: Dataset $D=(X_{1},\ldots,X_{m})$; parameters
$\varepsilon,\delta,\alpha,\beta\in(0,1)$

1:$\triangleright$ Some parameter settings.

2:$\phi\leftarrow 1$, $t\leftarrow 3/2$, $r\leftarrow 1$

3:$\eta\leftarrow\frac{\alpha}{C_{1}(d+\sqrt{d\ln(4/\beta)})}$
$\triangleright$ $C_{1}$ from Proposition 6.2

4:$\gamma\leftarrow\min\left\\{\sqrt{\frac{\varepsilon}{2d(d+1/\eta^{2})}},\frac{\varepsilon}{8d\sqrt{\ln(2/\delta)}},\frac{\varepsilon}{8\ln(2/\delta)},\frac{\varepsilon\eta}{12\sqrt{d}\sqrt{\ln(2/\delta)}}\right\\}$.

5:$k\leftarrow\max\left\\{\frac{400(r+\phi)}{\gamma},\frac{20}{\varepsilon}\ln\left(1+\frac{e^{\varepsilon}-1}{2\delta}\right)\right\\}$

6:function $\mathcal{A}$($X_{1},\ldots,X_{s},\alpha,\beta$)

7: $\widehat{\Sigma}=$
LearnCovariance$(\alpha,\beta/k,X_{1},\ldots,X_{s})$$\triangleright$ Algorithm
from [11]

8: Return $\widehat{\Sigma}$.

9:end function

10:function $\mathcal{B}$($\widehat{\Sigma}$)

11: Let $G$ be $d\times d$ matrix with independent $\mathcal{N}(0,1)$ entries.

12: Return $\widehat{\Sigma}^{1/2}(I+\eta G)(I+\eta
G)^{\top}\widehat{\Sigma}^{1/2}$.

13:end function

######  Theorem 6.3.

Applying 3.6 with parameters specified in Algorithm 4 gives a mechanism
$\mathcal{M}$ satisfying the following. For all $m,d\in\mathbb{N}$,
$\varepsilon,\delta,\alpha,\beta\in(0,1)$, and dataset
$D=(X_{1},\ldots,X_{m})$:

  1. 1.

Given $\varepsilon,\delta,\beta,\alpha,D$ as input, $\mathcal{M}$ is
$(2\varepsilon,4e^{\varepsilon}\delta)$-DP and runs in $poly(m,d)$ time.

  2. 2.

If $D$ is an $\alpha$-corrupted i.i.d. sample of size $m$ from
$\mathcal{N}(0,\Sigma)$, for $\Sigma\succ 0$, $\alpha$ is sufficiently small,
and
$m=\widetilde{\Omega}\left(\frac{d^{3.5}\ln^{5.5}(1/\beta)\ln(1/\delta)}{\varepsilon\alpha^{3}}\right)$,
then with probability at least $1-\beta$, $\mathcal{M}$ outputs a matrix
$\widehat{\Sigma}$ such that
$\|\widehat{\Sigma}^{-1/2}\Sigma\widehat{\Sigma}^{-1/2}-I_{d}\|_{F}=O(\alpha\ln(1/\alpha))$.

###### Proof.

The privacy claim follows from the exact same argument as in the proof of
Theorem 5.4 (still we have $400(r+\phi)/k\leq\gamma$). Also, the time
complexity of the algorithm remains polynomial since the non-private robust
covariance learning method runs in polynomial time.

We now prove utility although the proof does not directly make use of the
utility statement in 3.6. For $i\in[k]$, let $\widehat{\Sigma}_{i}$ be the
candidate covariance matrix computed by $\mathcal{A}$ with inputs
$\\{X_{\ell}\\}_{\ell=(i-1)s+1}^{is}$. Provided that
$s=\widetilde{\Omega}\left(\frac{d^{2}\log^{5}(k/\beta)}{\alpha^{2}}\right)$,
we have that $\|\Sigma^{-1/2}\widehat{\Sigma}_{i}\Sigma^{-1/2}-I_{d}\|_{F}\leq
C_{4}\alpha\ln(1/\alpha)<1/9$ for $i\in[k]$ with probability $1-\beta/2$. Note
that a sufficient condition for this is that $m\geq k\cdot
s=\widetilde{\Omega}\left(k\cdot\frac{d^{2}\log^{5}(k/\beta)}{\alpha^{2}}\right)=\widetilde{\Omega}\left(\frac{d^{3.5}\ln(1/\delta)\ln(1/\beta)^{5.5}}{\varepsilon\alpha^{3}}\right)$.
Next, using the trivial fact that Frobenius norm is an upper bound on spectral
norm, we also have
$\|\Sigma^{-1/2}\widehat{\Sigma}_{i}\Sigma^{-1/2}-I_{d}\|<1/9$. By A.6, we
also have
$\|\widehat{\Sigma}_{i}^{-1/2}\Sigma^{1/2}\widehat{\Sigma}_{i}^{-1/2}-I_{d}\|<4/9$.
Thus $\operatorname{dist}(\Sigma,\widehat{\Sigma}_{i})<4/9$, where
$\operatorname{dist}$ is as defined in Eq. (5.1). Now, using the fact that
$\operatorname{dist}(\cdot,\cdot)$ satisfies a $3/2$-approximate triangle
inequality, we have
$\operatorname{dist}(\widehat{\Sigma}_{i},\widehat{\Sigma}_{j})<2/3$ for all
$i,j\in[k]$. In particular, since $k$ is sufficiently large, the algorithm
does not fail with probability $1$ (conditioned on the event that
$\operatorname{dist}(\widehat{\Sigma}_{j},\widehat{\Sigma}_{j})<2/3$ for all
$i,j\in[k]$.

Now, let $\widehat{\Sigma}_{\mu}$ be the weighted average computed in Line 9
of Algorithm 1. Note that if
$\operatorname{dist}(\widehat{\Sigma}_{i},\widehat{\Sigma}_{j})<2/3$ for all
$i,j\in[k]$ then
$\widehat{\Sigma}_{\mu}=\frac{1}{k}\sum_{i\in[k]}\widehat{\Sigma}_{i}$. It is
easy to check that since
$\|\Sigma^{-1/2}\widehat{\Sigma}_{i}\Sigma^{-1/2}-I_{d}\|_{F}=O(\alpha\ln(1/\alpha))$,
we also have
$\|\Sigma^{-1/2}\widehat{\Sigma}_{\mu}\Sigma^{-1/2}-I_{d}\|_{F}=O(\alpha\ln(1/\alpha))$.

Next, note that by our choice of $\eta$ in Algorithm 4, Proposition 6.2 shows
that $\mathcal{B}$ is $(\alpha/\sqrt{d},\beta/2)$-concentrated. In other
words, if $\widehat{\Sigma}=\mathcal{B}(\widehat{\Sigma}_{\mu})$ then
$\|\widehat{\Sigma}^{-1/2}\widehat{\Sigma}_{\mu}\widehat{\Sigma}^{-1/2}-I_{d}\|=O(\alpha/\sqrt{d})$.
In particular,
$\|\widehat{\Sigma}^{-1/2}\widehat{\Sigma}_{\mu}\widehat{\Sigma}^{-1/2}-I_{d}\|_{F}=O(\alpha)$.

We are now ready to complete the proof. We have that
$d_{\mathrm{TV}}\left({\mathcal{N}(0,\Sigma)},{\mathcal{N}(0,\widehat{\Sigma}_{\mu})}\right)\leq
O(\alpha\ln(1/\alpha))$ and
$d_{\mathrm{TV}}\left({\mathcal{N}(0,\widehat{\Sigma}_{\mu})},{\mathcal{N}(0,\widehat{\Sigma})}\right)\leq
O(\alpha)$ by 2.1 and 2.2. By the triangle inequality, we conclude that
$d_{\mathrm{TV}}\left({\mathcal{N}(0,\Sigma)},{\mathcal{N}(0,\widehat{\Sigma})}\right)\leq
O(\alpha\ln(1/\alpha))$. ∎

##  7 Efficient and robust mean estimation

As another application of our general framework, we show how to efficiently
and robustly learn the mean of a Gaussian distribution with identity
covariance matrix. The following proposition indicates that every normed
vector space induces a convex semimetric space, which means the Euclidean
space is a convex semimetric space too.

######  Proposition 7.1.

Every normed vector space equipped with its induced metric forms a convex
semimetric space. In particular, it satisfies triangle inequality ($r=\infty$,
$t=1$) and $0$-locality.

###### Proof.

The first four properties of Definition 3.1 hold for norm-induced metrics by
standard arguments. Locality holds because

$\displaystyle\left\|\sum_{i\in[k]}\alpha_{i}Y_{i}-\sum_{i\in[k]}\alpha_{i}^{\prime}Y_{i}\right\|=$
$\displaystyle\left\|\sum_{i\in[k]}(\alpha_{i}-\alpha_{i}^{\prime})(Y_{i}-Y_{1})\right\|\leq$
$\displaystyle\sum_{i\in[k]}|\alpha_{i}-\alpha_{i}^{\prime}|\left\|(Y_{i}-Y_{1})\right\|\leq\sum_{i\in[k]}|\alpha_{i}-\alpha_{i}^{\prime}|\max_{i,j}\|Y_{i}-Y_{j}\|,$

where in the first equality we used that
$\sum_{i\in[k]}\alpha_{i}=\sum_{i\in[k]}\alpha_{i}^{\prime}=1$. ∎

The Gaussian mechanism is known to be a private mechanism for functions with
low $\ell_{2}$-sensitivity. The following theorem presents this observation in
the terminology of the masking mechanisms.

######  Lemma 7.2.

For $\mu\in\mathbb{R}^{d}$, define $\mathcal{B}(\mu)=\mu+\eta X$ where
$X\sim\mathcal{N}\left(0,I\right)$ and
$\eta=\frac{2\gamma\ln(1.25/\delta)}{\varepsilon}$. Then $\mathcal{B}$ is a
$(\gamma,\varepsilon,\delta)$-masking mechanism. Moreover, for all $\beta>0$,
$\mathcal{B}$ is $(\eta\sqrt{2d\ln(2/\beta)},\beta/2)$ concentrated.

###### Proof.

The privacy result is standard and is identical to that of the Gaussian
mechanism; see Theorem A.1. in [20] for details. The concentration result is
also standard, and follows directly from the Hoeffding’s inequality (see [45,
Proposition 2.5]). ∎

Now we have the ingredients to reduce the problem of privately learning the
mean to the non-private one. The following result is about the polynomial
learnability of the mean in the robust (non-private) setting and is a direct
consequence of Theorem 1.1 in [15].

######  Theorem 7.3 ([15]).

There an absolute constant $C_{3}$ and an algorithm that receives
$0<\alpha<1/2$, $\beta>0$, and $X_{1},\ldots,X_{m}$ as input, runs in
$poly(m,d)$, and outputs $\widehat{\mu}$ with the following guarantee: if
$X_{1},\ldots,X_{m}$ are $\alpha$-corrupted samples from
$\mathcal{N}(\mu,\Sigma)$ for some $(1/2)I_{d}\preceq\Sigma\preceq 2I_{d}$,
and if $m={\Omega}\left(\frac{d+\log 1/\beta}{\alpha^{2}}\right)$, then
$\left\|\mu-\widehat{\mu}\right\|_{2}\leq C_{3}\alpha\sqrt{\log(1/\alpha)}$.

Finally, we are ready to apply the 3.6. Algorithm 5 summarizes the parameters
that we use for this reduction.

Algorithm 5 A private and robust algorithm for learning the mean

Input: Dataset $D=(X_{1},\ldots,X_{m})$; parameters
$\varepsilon,\delta,\alpha,\beta\in(0,1)$

1:$\triangleright$ Some parameter settings.

2:$t\leftarrow 1,\phi\leftarrow 0$ $\triangleright$ See Proposition 7.1.

3:$r\leftarrow C_{3}\alpha\log(1/\alpha)$ $\triangleright$ $C_{3}$ is from
Theorem 7.3.

4:$\eta\leftarrow\frac{\alpha}{\sqrt{2d\ln(2/\beta)}}$

5:$\gamma\leftarrow\frac{\eta\varepsilon}{2\ln(1.25/\delta)}$

6:$k\leftarrow\max\left\\{\frac{400r}{\gamma},\frac{20}{\varepsilon}\ln\left(1+\frac{e^{\varepsilon}-1}{2\delta}\right)\right\\}$

7:function $\mathcal{A}$($X_{1},\ldots,X_{s},\alpha,\beta$)

8: $\widehat{\mu}=$
QUEScoreFilter$(\alpha,\beta/k,X_{1},\ldots,X_{s})$$\triangleright$ Algorithm
from [15]

9: Return $\widehat{\mu}$.

10:end function

11:function $\mathcal{B}$($\widehat{\mu}$)

12: Let $X\sim\mathcal{N}(0,I_{d})$

13: Return $\widehat{\mu}+\eta X$.

14:end function

######  Theorem 7.4.

Applying 3.6 with parameters specified in Algorithm 5 gives a mechanism
$\mathcal{M}$ satisfying the following. For all $m,d\in\mathbb{N}$,
$\varepsilon,\delta,\alpha,\beta\in(0,1)$, and dataset
$D=(X_{1},\ldots,X_{m})$:

  1. 1.

Given $\varepsilon,\delta,\beta,\alpha,D$ as input, $\mathcal{M}$ is
$(2\varepsilon,4e^{\varepsilon}\delta)$-DP and runs in $poly(m,d)$ time.

  2. 2.

If $D$ is an $\alpha$-corrupted i.i.d. sample of size $m$ from
$\mathcal{N}(\mu,\Sigma)$, for $(1/2)I\preceq\Sigma\preceq 2I$ ,
$\alpha\in(0,0.5)$, and
$m=\widetilde{\Omega}\left(\frac{d^{3/2}\ln(1/\delta)}{\varepsilon\alpha^{2}}\right)$,
then with probability at least $1-\beta$, $\mathcal{M}$ outputs a vector
$\widetilde{\mu}$ such that
$\|\widetilde{\mu}-\mu\|_{2}=O(\alpha\sqrt{\ln(1/\alpha)})$.

###### Proof.

Proposition 7.1 indicates that $(\mathbb{R}^{d},\|.\|_{2})$ form a convex
semimetric space with non-restricted triangle inequality and $0$-locality. We
take $r=C_{3}\alpha\sqrt{\log(1/\alpha)}<\infty$. By 7.2, $\mathcal{B}$ is
$(\gamma,\varepsilon,\delta)$-masking mechanism with $\gamma,\eta$ as defined
in Algorithm 5. Our choice of $k$ in Algorithm 5 means that $400r/k\leq\gamma$
so $\mathcal{B}$ is also a $(400r/k,\varepsilon,\delta)$-masking mechanism.
Therefore, the privacy of the algorithm follows from 3.6. Also, the time
complexity of the algorithm remains polynomial since the non-private robust
mean estimation runs in polynomial time.

7.2 shows that $\mathcal{B}$ is
$(\eta\sqrt{2d\ln(2/\beta)},\beta/2)$-concentrated. By our choice of $\eta$ in
Algorithm 5, we conclude that $\mathcal{B}$ is
$(\alpha,\beta/2)$-concentrated.

Let $\mu_{i}=\mathcal{A}(\\{X_{\ell}\\}^{is}_{\ell=(i-1)s+1},\alpha,\beta)$ be
the outputs of the non-private mean estimation method. Since
$m=\widetilde{\Omega}\left(k\cdot\frac{d+\log(k/\beta)}{\alpha^{2}}\right)$,
we conclude that $\|\mu_{i}-\mu\|_{2}=O(\alpha\sqrt{\log(1/\alpha)})$ for all
$i\in[k]$ with probability $1-\beta/2$ by Theorem 7.3 and a union bound over
$k$.

Finally, setting $Y_{i}=\mu_{i}$, $Y^{*}=\mu$, and
$\widehat{Y}=\widetilde{\mu}$ for Algorithm 1, with the modules specified in
Algorithm 4. Then, taking
$\alpha_{1}=\alpha_{2}=O(\alpha\sqrt{\log(1/\alpha)})$ in 3.6 gives that
$\operatorname{dist}(\widetilde{\mu},\mu)=O(\alpha\sqrt{\log(1/\alpha)})$ with
probability $1-\beta$. Moreover, it is enough for the number of samples to be
$m=\widetilde{\Omega}\left(k\cdot\frac{d+\log(k/\beta)}{\alpha^{2}}\right)=\widetilde{\Omega}\left(\frac{d^{3/2}\ln(1/\delta)}{\varepsilon\alpha^{2}}\right)$,
which concludes the proof. ∎

Finally, we can combine the mean estimation approach of Theorem 7.4 with the
covariance estimation of Theorem 6.3 to have a complete private algorithm for
robust learning of Gaussians. The sample complexity bottleneck will be that of
covariance estimation. The proof of the following corollary is very similar to
that of Subsection 5.1 and is omitted for brevity.

inline,color=red!30]Chris: Check, below should probably be $d^{3.5}$.

######  Corollary 7.5.

There exist an $(O(\varepsilon),O(e^{\varepsilon}\delta))$-DP algorithm that
receives $0<\alpha<1/2$, $\beta>0$, and $X_{1},\ldots,X_{m}$ as input, and
outputs $\widehat{\Sigma},\widehat{\mu}$ with the following guarantee: if
$X_{1},\ldots,X_{m}$ are $\alpha$-corrupted samples from
$\mathcal{N}(\mu,\Sigma)$ for
$\Sigma\in\mathbb{S}_{>0}^{d},\mu\in\mathbb{R}^{d}$ then
$d_{TV}(\mathcal{N}(\mu,\Sigma),\mathcal{N}(\widehat{\mu},\widehat{\Sigma}))=O(\alpha\log(1/\alpha))$
as long as
$m={\Omega}\left(\frac{d^{4.5}\ln(1/\delta)\ln^{5}(d/\varepsilon\alpha\beta)}{\varepsilon\alpha^{3}}\right)$.

######  Remark 7.6.

Unlike Theorem 5.7 for the non-robust case, here the utility guarantee is
conditioned on the assumption that $\Sigma\succ 0$. In order to address this,
one needs to come up with either a variant of the robust covariance estimation
approach of [11] which works for singular covariance matrices, or an efficient
and robust subspace recovery method which can be a hard task [23].

######  Remark 7.7.

The sample complexities of Theorem 6.3 (and therefore 7.5) as well as Theorem
7.4 are not optimal. It seems possible to improve these sample complexities by
a more careful analysis. We leave this as a future research direction.

## References

  * [1] Ishaq Aden-Ali, Hassan Ashtiani, and Gautam Kamath.  On the sample complexity of privately learning unbounded high-dimensional gaussians.  In Algorithmic Learning Theory, pages 185–216. PMLR, 2021.
  * [2] Ishaq Aden-Ali, Hassan Ashtiani, and Christopher Liaw.  Privately learning mixtures of axis-aligned gaussians.  Advances in Neural and Information Processing Systems, 2021.
  * [3] Hassan Ashtiani, Shai Ben-David, Nicholas JA Harvey, Christopher Liaw, Abbas Mehrabian, and Yaniv Plan.  Near-optimal sample complexity bounds for robust learning of gaussian mixtures via compression schemes.  Journal of the ACM (JACM), 67(6):1–42, 2020.
  * [4] Hassan Ashtiani and Christopher Liaw.  Private and polynomial time algorithms for learning gaussians and beyond: Extended abstract.  Conference on Learning Theory, 2022.
  * [5] Sourav Biswas, Yihe Dong, Gautam Kamath, and Jonathan Ullman.  Coinpress: Practical private mean and covariance estimation.  Advances in Neural and Information Processing Systems, 2020.
  * [6] Gavin Brown, Marco Gaboardi, Adam Smith, Jonathan Ullman, and Lydia Zakynthinou.  Covariance-aware private mean estimation without private covariance estimation.  arXiv preprint arXiv:2106.13329, 2021.
  * [7] Mark Bun, Gautam Kamath, Thomas Steinke, and Steven Z Wu.  Private hypothesis selection.  Advances in Neural Information Processing Systems, 32:156–167, 2019\.
  * [8] Mark Bun, Kobbi Nissim, and Uri Stemmer.  Simultaneous private learning of multiple concepts.  Journal of Machine Learning Research, 20(94):1–34, 2019.
  * [9] Mark Bun and Thomas Steinke.  Concentrated differential privacy: Simplifications, extensions, and lower bounds.  In Theory of Cryptography Conference, pages 635–658. Springer, 2016\.
  * [10] Kamalika Chaudhuri, Anand D Sarwate, and Kaushik Sinha.  A near-optimal algorithm for differentially-private principal components.  Journal of Machine Learning Research, 14, 2013.
  * [11] Ilias Diakonikolas, Gautam Kamath, Daniel Kane, Jerry Li, Ankur Moitra, and Alistair Stewart.  Robust estimators in high-dimensions without the computational intractability.  SIAM Journal on Computing, 48(2):742–864, 2019.
  * [12] Ilias Diakonikolas, Gautam Kamath, Daniel M Kane, Jerry Li, Ankur Moitra, and Alistair Stewart.  Robust estimators in high dimensions without the computational intractability.  In 2016 IEEE 57th Annual Symposium on Foundations of Computer Science (FOCS), pages 655–664. IEEE, 2016.
  * [13] Wei Dong, Yuting Liang, and Ke Yi.  Differentially private covariance revisited.  arXiv preprint arXiv:2205.14324, 2022.
  * [14] Wei Dong and Ke Yi.  Universal private estimators.  arXiv preprint arXiv:2111.02598, 2021.
  * [15] Yihe Dong, Samuel Hopkins, and Jerry Li.  Quantum entropy scoring for fast robust mean estimation and improved outlier detection.  Advances in Neural Information Processing Systems, 32:6067–6077, 2019.
  * [16] Wenxin Du, Canyon Foot, Monica Moniot, Andrew Bray, and Adam Groce.  Differentially private confidence intervals.  arXiv preprint arXiv:2001.02285, 2020.
  * [17] Cynthia Dwork, Krishnaram Kenthapadi, Frank McSherry, Ilya Mironov, and Moni Naor.  Our data, ourselves: Privacy via distributed noise generation.  In Annual International Conference on the Theory and Applications of Cryptographic Techniques, pages 486–503. Springer, 2006.
  * [18] Cynthia Dwork and Jing Lei.  Differential privacy and robust statistics.  In Proceedings of the forty-first annual ACM symposium on Theory of computing, pages 371–380, 2009.
  * [19] Cynthia Dwork, Frank McSherry, Kobbi Nissim, and Adam Smith.  Calibrating noise to sensitivity in private data analysis.  Journal of Privacy and Confidentiality, 7(3):17–51, 2017.
  * [20] Cynthia Dwork and Aaron Roth.  The algorithmic foundations of differential privacy.  Found. Trends Theor. Comput. Sci., 9(3-4):211–407, 2014.
  * [21] Cynthia Dwork and Guy N Rothblum.  Concentrated differential privacy.  arXiv preprint arXiv:1603.01887, 2016.
  * [22] Quan Geng, Wei Ding, Ruiqi Guo, and Sanjiv Kumar.  Tight analysis of privacy and utility tradeoff in approximate differential privacy.  Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, 108:89–99, 2020.
  * [23] Moritz Hardt and Ankur Moitra.  Algorithms and hardness for robust subspace recovery.  In Conference on Learning Theory, pages 354–375. PMLR, 2013.
  * [24] Gautam Kamath, Jerry Li, Vikrant Singhal, and Jonathan Ullman.  Privately learning high-dimensional distributions.  In Conference on Learning Theory, pages 1853–1902. PMLR, 2019.
  * [25] Gautam Kamath, Argyris Mouzakis, Vikrant Singhal, Thomas Steinke, and Jonathan Ullman.  A private and computationally-efficient estimator for unbounded gaussians.  arXiv preprint arXiv:2111.04609, 2021.
  * [26] Gautam Kamath, Or Sheffet, Vikrant Singhal, and Jonathan Ullman.  Differentially private algorithms for learning mixtures of separated gaussians.  In 2020 Information Theory and Applications Workshop (ITA), pages 1–62. IEEE, 2020.
  * [27] Gautam Kamath, Vikrant Singhal, and Jonathan Ullman.  Private mean estimation of heavy-tailed distributions.  In Conference on Learning Theory, pages 2204–2235. PMLR, 2020.
  * [28] Gautam Kamath and Jonathan Ullman.  A primer on private statistics.  arXiv preprint arXiv:2005.00010, 2020.
  * [29] Vishesh Karwa and Salil Vadhan.  Finite sample differentially private confidence intervals.  In 9th Innovations in Theoretical Computer Science Conference (ITCS 2018). Schloss Dagstuhl-Leibniz-Zentrum fuer Informatik, 2018.
  * [30] Pravesh K Kothari, Pasin Manurangsi, and Ameya Velingker.  Private robust estimation by stabilizing convex relaxations.  arXiv preprint arXiv:2112.03548, 2021.
  * [31] Solomon Kullback.  Information Theory and Statistics.  Dover, 1968.
  * [32] Kevin A Lai, Anup B Rao, and Santosh Vempala.  Agnostic estimation of mean and covariance.  In 2016 IEEE 57th Annual Symposium on Foundations of Computer Science (FOCS), pages 665–674. IEEE, 2016.
  * [33] B. Laurent and P. Massart.  Adaptive estimation of a quadratic functional by model selection.  The Annals of Statistics, 28(5):1302–1338, 2000.
  * [34] Xiyang Liu, Weihao Kong, Prateek Jain, and Sewoong Oh.  Dp-pca: Statistically optimal and differentially private pca.  arXiv preprint arXiv:2205.13709, 2022.
  * [35] Xiyang Liu, Weihao Kong, Sham Kakade, and Sewoong Oh.  Robust and differentially private mean estimation.  arXiv preprint arXiv:2102.09159, 2021.
  * [36] Kobbi Nissim, Sofya Raskhodnikova, and Adam Smith.  Smooth sensitivity and sampling in private data analysis.  In Proceedings of the thirty-ninth annual ACM symposium on Theory of computing, pages 75–84, 2007.
  * [37] Kobbi Nissim, Uri Stemmer, and Salil Vadhan.  Locating a small cluster privately.  In Proceedings of the 35th ACM SIGMOD-SIGACT-SIGAI Symposium on Principles of Database Systems, pages 413–427, 2016.
  * [38] Carl Edward Rasmussen and Christopher K. I. Williams.  Gaussian Processes for Machine Learning.  The MIT Press, 2006.
  * [39] Vikrant Singhal and Thomas Steinke.  Privately learning subspaces.  arXiv preprint arXiv:2106.00001, 2021.
  * [40] Adam Smith.  Privacy-preserving statistical estimation with optimal convergence rates.  In Proceedings of the forty-third annual ACM symposium on Theory of computing, pages 813–822, 2011.
  * [41] Eliad Tsfadia, Edith Cohen, Haim Kaplan, Yishay Mansour, and Uri Stemmer.  Friendlycore: Practical differentially private aggregation.  arXiv preprint arXiv:2110.10132, 2021.
  * [42] Alexandre B Tsybakov.  Introduction to Nonparametric Estimation. Springer, 2009.
  * [43] Salil Vadhan.  The complexity of differential privacy.  In Tutorials on the Foundations of Cryptography, pages 347–450. Springer, 2017.
  * [44] Roman Vershynin.  High-Dimensional Probability: An Introduction with Applications in Data Science.  Cambridge University Press, 2018.
  * [45] Martin J. Wainwright.  High-Dimensional Statistics: A Non-Asymptotic Viewpoint.  Cambridge University Press, 2019.
  * [46] Lu Wei, Anand D Sarwate, Jukka Corander, Alfred Hero, and Vahid Tarokh.  Analysis of a privacy-preserving pca algorithm using random matrix theory.  In 2016 IEEE Global Conference on Signal and Information Processing (GlobalSIP), pages 1335–1339. IEEE, 2016.

##  Appendix A Standard Facts

######  Claim A.1.

If $x\leq 1/2$ then $\frac{1+x}{1-x}\leq 1+4x$.

###### Proof.

Let $f(x)=1+4x-\frac{1+x}{1-x}$. Then
$f^{\prime\prime}(x)=-\frac{2(1+x)}{(1-x)^{3}}-\frac{2}{(1-x)^{2}}$. So $f$ is
concave on the interval $[0,1)$. Since $f(0)=f(1/2)=0$, we have that $f(x)\geq
0$ for $x\in[0,1/2]$ as desired. ∎

######  Claim A.2.

For all $x\in\mathbb{R}$, $\frac{1-x}{1+x}\geq 1-2x$.

###### Proof.

We have $\frac{1-x}{1+x}\geq(1-x)^{2}=1-2x+x^{2}\geq 1-2x$. ∎

######  Lemma A.3 ([33, Lemma 1]).

Let $g_{1},\ldots,g_{k}$ be i.i.d. $\mathcal{N}(0,1)$ random variables. Then

$\mathbb{P}\left[\sum_{i=1}^{k}g_{i}^{2}-k\geq 2\sqrt{kt}+2t\right]\leq
e^{-t}.$

######  Lemma A.4 ([44, Exercise 4.7.3]).

There is an absolute constant $C>0$ such that the following holds. Let
$\Sigma\succ 0$ and $\beta>0$. Let
$X_{1},\ldots,X_{m}\sim\mathcal{N}(0,\Sigma)$ and
$\widehat{\Sigma}=\frac{1}{m}\sum_{i=1}^{m}X_{i}X_{i}^{\top}$. Then
$\|\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2}-I_{d}\|\leq
C\left(\sqrt{\frac{d+\ln(2/\beta)}{m}}+\frac{d+\ln(2/\beta)}{m}\right)$ with
probability $1-\beta$.

######  Theorem A.5 ([44, Theorem 4.4.5]).

There is a universal constant $C>0$ such that the following holds. Let $G$ be
a $d\times d$ random matrix where each entry is an independent
$\mathcal{N}(0,1)$ random variable. Then $\|G\|\leq
C(\sqrt{d}+\sqrt{\ln(2/\beta)})$ with probability $1-\beta$.

######  Lemma A.6.

Let $A,B$ be $d\times d$ matrices and suppose that
$\|A^{-1/2}BA^{-1/2}-I_{d}\|\leq\gamma\leq 1/2$. Then
$\|B^{-1/2}AB^{-1/2}-I_{d}\|\leq 4\gamma$.

###### Proof.

The assumption that $\|A^{-1/2}BA^{-1/2}-I_{d}\|\leq\gamma$ means that
$(1-\gamma)A\preceq B\preceq(1+\gamma)A$. Hence, we also have that
$\frac{1}{1+\gamma}B\preceq A\preceq\frac{1}{1-\gamma}B$. The asserted
inequality is then proved using the bounds in Claim A.7. ∎

######  Claim A.7.

Suppose that $\gamma\in[0,1/2]$. Then $\frac{1}{1+\gamma}\geq 1-4\gamma$ and
$\frac{1}{1-\gamma}\leq 1+4\gamma$.

###### Proof.

To prove the first inequality, let $f(\gamma)=1-4\gamma-\frac{1}{1+\gamma}$.
We need to show that $f(\gamma)\geq 0$ for all $\gamma\in[0,1/2]$.
Differentiating, we have $f^{\prime}(\gamma)=4-\frac{1}{1+x^{2}}$ so $f$ is
increasing. Since $f(0)=0$, the first inequality follows.

For the second inequality, let $f(\gamma)=1+4\gamma-\frac{1}{1-\gamma}$. The
derivative $f^{\prime}(\gamma)=4-\frac{1}{(1-\gamma)^{2}}\geq 0$ for
$\gamma\in[0,1/2]$. So $f$ is increasing on $[0,1/2]$ and $f(0)=0$ so the
second inequality is proved. ∎

##  Appendix B Missing Proofs From Section 2

###  B.1 Proof of 2.1

###### Proof of 2.1.

The first assertion is a well-known formula for the KL divergence between two
normal distributions (see [31, §9.1] or [38, Eq. A.23]).

Next, let $\lambda_{1},\ldots,\lambda_{d}$ be the eigenvalues of
$\Sigma_{2}^{-1}\Sigma_{1}$. Then

$\displaystyle\operatorname{tr}(\Sigma_{2}^{-1}\Sigma_{1}-I)-\ln\det(\Sigma_{2}^{-1}\Sigma_{1})=\sum_{i=1}^{d}\lambda_{i}-1-\ln\lambda_{i}\leq\sum_{i=1}^{d}(\lambda_{i}-1)^{2}=\|\Sigma_{2}^{-1}\Sigma_{1}-I\|_{F}^{2},$

where the inequality follows from Claim B.1. The lemma now follows because the
two matrices $\Sigma_{2}^{-1}\Sigma_{1}$ and
$\Sigma_{2}^{-1/2}\Sigma_{1}\Sigma_{2}^{-1/2}$ have the same spectrum. ∎

######  Claim B.1.

Suppose that $x\geq 1/2$. Then $x-1-\ln(x)\leq(x-1)^{2}$.

###### Proof.

Let $f(x)=(x-1)^{2}-(x-1)+\ln(x)$. Then
$f^{\prime}(x)=2(x-1)-1+\frac{1}{x}=\frac{2x^{2}-3x+1}{x}=\frac{(2x-1)(x-1)}{x}$.
So $f^{\prime}(x)\leq 0$ on $[1/2,1]$ and $f^{\prime}(x)\geq 0$ on
$[1,\infty)$. Hence, on the interval $[1/2,\infty)$, $f(x)$ is minimized at
$x=1$ and $f(x)\geq f(1)=0$ for all $x\in[1/2,\infty)$. ∎

######  Claim B.2.

Suppose $x\leq 1/2$. Then $\sqrt{1-x}\geq 1-x/2-x^{2}/2$.

###### Proof.

Let $f(x)=\sqrt{1-x}-1+x/2+x^{2}/2$. Then
$f^{\prime}(x)=x-\frac{1}{2\sqrt{1-x}}+\frac{1}{2}$ and
$f^{\prime\prime}(x)=1-\frac{1}{4(1-x)^{3/2}}$. Note that $f^{\prime}(0)=0$ so
$x=0$ is a critical point of $f$. Moreover, $f^{\prime\prime}(x)\geq 0$ for
all $x\leq 1/2$ so $f$ is convex on $(-\infty,1/2)$. Hence, the minimizer of
$f$ on the interval $(-\infty,1/2]$ is at $x=0$ and $f(x)\geq f(0)=0$ for all
$x\in(-\infty,1/2]$. ∎

###  B.2 Proof of 2.4

######  Claim B.3 (see [45, Example 2.8]).

Let $X\sim\mathcal{N}(0,1)$ and $|u|\leq 1/4$. Then
$\mathbb{E}\left[e^{u(X^{2}-1)}\right]\leq e^{2u^{2}}$.

We first prove 2.4 for the special case where the matrix is diagonal. The
general case follows easily from rotation invariance.

######  Lemma B.4.

Let $\Lambda$ be a $d\times d$ diagonal matrix and
$g\sim\mathcal{N}(0,I_{d})$. Then for all $x>0$,

$\mathbb{P}\left[g^{\top}\Lambda
g\geq\operatorname{tr}(\Lambda)+4\|\Lambda\|_{F}\cdot\sqrt{x}+4\|\Lambda\|\cdot
x\right]\leq e^{-x}.$

###### Proof.

Note that $g^{\top}\Lambda g=\sum_{i=1}^{d}\lambda_{i}(g_{i}^{2}-1)$ where
$\lambda_{1},\ldots,\lambda_{d}$ are the diagonal entries of $\Lambda$. From
Claim B.3, we have $\mathbb{E}\left[e^{u\lambda_{i}(g_{i}^{2}-1)}\right]\leq
e^{2u^{2}\lambda_{i}^{2}}$ provided $|u|\leq 1/4\|\Lambda\|$. Hence,
$\mathbb{E}\left[e^{u(g^{\top}\Lambda
g-\operatorname{tr}(\Lambda))}\right]=\prod_{i=1}^{d}\mathbb{E}\left[e^{u\lambda_{i}(g_{i}^{2}-1)}\right]\leq
e^{2u^{2}\|\Lambda\|_{F}^{2}}$. Applying Markov’s Inequality, we thus have,
for $0\leq u\leq 1/4\|\Lambda\|$,

$\displaystyle\mathbb{P}\left[g^{\top}\Lambda g-\operatorname{tr}(\Lambda)\geq
4\|\Lambda\|_{F}\cdot\sqrt{x}+4\|\Lambda\|\cdot x\right]$ $\displaystyle\leq
e^{-u(4\|\Lambda\|_{F}\cdot\sqrt{x}+4\|\Lambda\|\cdot
x)+2u^{2}\|\Lambda\|_{F}^{2}}.$

It is a straightforward calculation to verify that choosing
$u=\min\left\\{\frac{1}{4\|\Lambda\|},\frac{\sqrt{x}}{\|\Lambda\|_{F}}\right\\}$
makes the RHS at most $e^{-x}$ which proves the lemma. ∎

###### Proof of 2.4.

Let $A=U\Lambda U^{\top}$ be the eigendecomposition of $A$. Then
$g^{\top}Ag=g^{\top}U\Lambda U^{\top}g\stackrel{{\scriptstyle
d}}{{=}}g^{\top}\Lambda g$ where the last equality in distribution is by
rotation invariance (B.5). The claim now follows from B.4 since $A$ and
$\Lambda$ have the same spectrum. ∎

######  Fact B.5.

Let $g\sim\mathcal{N}(0,I_{d})$ and $U$ be a $d\times d$ orthogonal matrix.
Then $Ug\sim\mathcal{N}(0,I_{d})$.

###  B.3 Proof of 2.10

###### Proof of 2.10.

Let $f_{1},f_{2}$ be the densities of $\mathcal{D}_{1},\mathcal{D}_{2}$,
respectively. Let $S\subseteq\mathbb{R}^{d}$ and let
$T=\\{y\in\mathbb{R}^{d}\,:\,f_{1}(y)\geq e^{\varepsilon}f_{2}(y)\\}$. The
assumption
$\mathbb{P}_{Y\sim\mathcal{D}_{1}}\left[\mathcal{L}_{{\mathcal{D}_{1}}\parallel{\mathcal{D}_{2}}}(Y)\geq\varepsilon\right]\leq\delta$
means that $\mathbb{P}_{Y\sim\mathcal{D}_{1}}\left[Y\in T\right]\leq\delta$.
Thus, we have

$\displaystyle\mathbb{P}_{Y\sim\mathcal{D}_{1}}\left[Y\in S\right]$
$\displaystyle=\int_{S}f_{1}(y)\,\mathrm{d}y$ $\displaystyle=\int_{S\setminus
T}f_{1}(y)\,\mathrm{d}y+\int_{S\cap T}f_{1}(y)\,\mathrm{d}y$
$\displaystyle\leq\int_{S\setminus
T}e^{\varepsilon}f_{2}(y)\,\mathrm{d}y+\delta$ $\displaystyle\leq
e^{\varepsilon}\mathbb{P}_{Y\sim\mathcal{D}_{2}}\left[S\right]+\delta.$

An identical argument shows that $\mathbb{P}_{Y\sim\mathcal{D}_{2}}\left[Y\in
S\right]\leq
e^{\varepsilon}\mathbb{P}_{Y\sim\mathcal{D}_{1}}\left[S\right]+\delta$, which
proves that $\mathcal{D}_{1},\mathcal{D}_{2}$ are
$(\varepsilon,\delta)$-indistinguishable. ∎

##  Appendix C Missing proofs from Section 3

###  C.1 Proof of 3.2

###### Proof.

The first two conditions in Definition 3.1 are trivial.

Triangle Inequality.

Next, we prove that $\operatorname{dist}$ satisfies a $(3/2)$-approximate
$1$-restricted triangle inequality. Let $A,B,C\in\mathcal{S}^{d}$ and let
$\rho_{ab}=\operatorname{dist}(A,B)$, $\rho_{bc}=\operatorname{dist}(B,C)$ and
$\rho_{ac}=\operatorname{dist}(A,C)$. We need to show that if
$\rho_{ab},\rho_{bc}\leq 1$ then
$\rho_{ac}\leq\frac{3}{2}(\rho_{ab}+\rho_{bc})$. Indeed, we have

$\displaystyle A$
$\displaystyle\preceq(1+\rho_{ab})B\preceq(1+\rho_{ab})(1+\rho_{bc})C=(1+\rho_{ab}+\rho_{bc}+\rho_{ab}\rho_{bc})C$
$\displaystyle\preceq(1+\rho_{ab}+\rho_{bc}+\rho_{ab}^{2}/2+\rho_{bc}^{2}/2)C\preceq(1+\frac{3}{2}(\rho_{ab}+\rho_{bc}))C,$

where the first two inequalities used that
$\operatorname{dist}(A,B)=\rho_{ab}$ and $\operatorname{dist}(B,C)=\rho_{bc}$,
respectively, the third inequality used that $xy\leq x^{2}/2+y^{2}/2$ for any
$x,y\in\mathbb{R}$ and the last inequality used that $\rho_{ab},\rho_{bc}\leq
1$. Similarly, we can show that

$A\succeq(1-\frac{3}{2}(\rho_{ab}+\rho_{bc}))C,$

which imply, with the previous quality, that
$\|C^{-1/2}AC^{-1/2}-I\|\leq\frac{3}{2}(\rho_{ab}+\rho_{bc})$. We can
similarly show that
$\|A^{-1/2}CA^{-1/2}-I\|\leq\frac{3}{2}(\rho_{ab}+\rho_{bc})$ so
$\rho_{ac}\leq\frac{3}{2}(\rho_{ab}+\rho_{bc})$, as desired.

Convexity.

For this, it would suffice to show that for all $\alpha\in[0,1]$,
$\operatorname{dist}(\alpha
A+(1-\alpha)B,C)\leq\alpha\operatorname{dist}(A,C)+(1-\alpha)\operatorname{dist}(B,C)$,
and the rest follows from induction over $k$. Let
$\rho_{a}=\operatorname{dist}(A,C),\rho_{b}=\operatorname{dist}(B,C)$. We need
to show the following four inequalities:

  1. 1.

$C\preceq(1+\alpha\rho_{a}+(1-\alpha)\rho_{b})(\alpha A+(1-\alpha)B)$;

  2. 2.

$C\succeq(1-\alpha\rho_{a}-(1-\alpha)\rho_{b})(\alpha A+(1-\alpha)B)$;

  3. 3.

$\alpha A+(1-\alpha)B\preceq(1+\alpha\rho_{a}+(1-\alpha)\rho_{b})C$; and

  4. 4.

$\alpha A+(1-\alpha)B\succeq(1-\alpha\rho_{a}-(1-\alpha)\rho_{b})C$.

The third and fourth items are easy to check. Indeed, for the third item, we
have $A\preceq(1+\rho_{a})C$ and $B\preceq(1+\rho_{b})C$ so

$\alpha
A+(1-\alpha)B\preceq\alpha(1+\rho_{a})C+(1-\alpha)(1+\rho_{b})C=(1+\alpha\rho_{a}+(1-\alpha)\rho_{b})C.$

The fourth item follows similarly.

It remains to prove the first two inequalities. Here, we prove only the first
inequality since the second follows by an analogous argument. Let

$K=\left(\frac{\alpha(1+\alpha\rho_{a}+(1-\alpha)\rho_{b})}{1+\rho_{a}}+\frac{(1-\alpha)(1+\alpha\rho_{a}+(1-\alpha)\rho_{b})}{1+\rho_{b}}\right)^{-1}$

and set

$w=K\frac{\alpha(1+\alpha\rho_{a}+(1-\alpha)\rho_{b})}{1+\rho_{a}}.$

Our choice of $K$ means that

$1-w=K\frac{(1-\alpha)(1+\alpha\rho_{a}+(1-\alpha)\rho_{b})}{1+\rho_{b}}.$

Note that $w\in[0,1]$. Since $C\preceq(1+\rho_{a})A$ and
$C\preceq(1+\rho_{b})B$, we have

$\displaystyle C$ $\displaystyle=wC+(1-w)C$ $\displaystyle\preceq
w(1+\rho_{a})C+(1-w)(1+\rho_{b})B$
$\displaystyle=K\left[\alpha(1+\alpha\rho_{a}+(1-\alpha)\rho_{b})A+(1-\alpha)(1+\alpha\rho_{a}+(1-\alpha)\rho_{b})B\right]$
$\displaystyle=K\left[(1+\alpha\rho_{a}+(1-\alpha)\rho_{b})(\alpha
A+(1-\alpha)B)\right].$

Claim C.1 shows that $K\leq 1$ so we have the desired inequality.

######  Claim C.1.

Let $\alpha\in[0,1]$ and let $\rho_{a},\rho_{b}\in\mathbb{R}_{\geq 0}$. Then

$\frac{\alpha(1+\alpha\rho_{a}+(1-\alpha)\rho_{b})}{1+\rho_{a}}+\frac{(1-\alpha)(1+\alpha\rho_{a}+(1-\alpha)\rho_{b})}{1+\rho_{b}}\geq
1.$

###### Proof.

Let

$f(\alpha)\coloneqq\frac{\alpha(1+\alpha\rho_{a}+(1-\alpha)\rho_{b})}{1+\rho_{a}}+\frac{(1-\alpha)(1+\alpha\rho_{a}+(1-\alpha)\rho_{b})}{1+\rho_{b}}.$

Observe that $f(0)=f(1)=1$. Thus, it suffices to prove that $f$ is concave.
Indeed, we have

$f^{\prime}(\alpha)=\frac{1+\alpha\rho_{a}+(1-\alpha)\rho_{b}}{1+\rho_{a}}+\frac{\alpha(\rho_{a}-\rho_{b})}{1+\rho_{a}}-\frac{1+\alpha\rho_{a}+(1-\alpha)\rho_{b}}{1+\rho_{b}}+\frac{(1-\alpha)(\rho_{a}-\rho_{b})}{1+\rho_{b}}.$

and

$f^{\prime\prime}(\alpha)=2\frac{\rho_{a}-\rho_{b}}{1+\rho_{a}}+2\frac{\rho_{b}-\rho_{a}}{1+\rho_{b}}=2(\rho_{a}-\rho_{b})\cdot\left(\frac{1}{1+\rho_{a}}-\frac{1}{1+\rho_{b}}\right).$

Finally, we check that $f^{\prime\prime}\leq 0$. Indeed, if
$\rho_{a}=\rho_{b}$ then $f^{\prime\prime}=0$ and if $\rho_{1}\neq\rho_{b}$
then $\rho_{a}-\rho_{b}$ and $1/(1+\rho_{a})-1/(1+\rho_{b})$ have opposite
signs so $f^{\prime\prime}<0$. ∎

Locality.

Let $\widehat{\Sigma}=\sum_{i\in[k]}w_{i}\widehat{\Sigma}_{i}$ and
$\widehat{\Sigma}^{\prime}=\sum_{i\in[k]}w_{i}^{\prime}\widehat{\Sigma}_{i}$.
Let
$\rho=\max_{i,j\in[k]}\operatorname{dist}(\widehat{\Sigma}_{i},\widehat{\Sigma}_{j})$.
Then

$\displaystyle\sum_{i\in[k]}w_{i}\widehat{\Sigma}_{i}$
$\displaystyle\preceq\sum_{i\in[k]}w_{i}^{\prime}\widehat{\Sigma}_{i}+\sum_{i\in[k]}|w_{i}-w_{i}^{\prime}|\widehat{\Sigma}_{i}$
$\displaystyle\preceq\widehat{\Sigma}^{\prime}+\sum_{i\in[k]}|w_{i}-w_{i}^{\prime}|(1+\rho)\widehat{\Sigma}^{\prime}$
$\displaystyle=\left(1+(1+\rho)\sum_{i\in[k]}|w_{i}-w_{i}^{\prime}|\right)\widehat{\Sigma}^{\prime}$

For the second inequality, we used
$\widehat{\Sigma}_{i}\preceq(1+\rho)\widehat{\Sigma}_{j}$ for all $i,j\in[k]$,
which imply $\widehat{\Sigma}_{i}\preceq(1+\rho)\widehat{\Sigma}^{\prime}$ for
all $i\in[k]$. ∎

###  C.2 Proof of 3.6

###### Proof.

We first prove the utility and then the privacy result.

#### Utility.

We first show that the algorithm will not return $\perp$ if the assumptions of
the utility guarantee hold. Based on the assumption,
$\operatorname{dist}(Y_{i},Y^{*})\leq\alpha_{1}/t$ for all $i\in[k]$. By the
$t$-approximate $r$-restricted triangle inequality we have
$\operatorname{dist}(Y_{i},Y_{j})\leq\alpha_{1}\leq r/t$ for all $i,j\in[k]$.
Therefore, for all $i\in[k]$, $q_{i}=1$. Finally,
$Q=1>0.8+\frac{2}{k\varepsilon}\ln{\left(1+\frac{e^{\varepsilon}-1}{2\delta}\right)}$
and therefore the algorithm does not fail in Line 7.

Let $W=\sum_{i\in[k]}{w_{i}}$. To show that $\widehat{Y}$ concentrates around
$Y^{*}$ we have

$\displaystyle\mathbb{P}\left[\operatorname{dist}(\hat{Y},Y^{*})>\alpha_{1}+\alpha_{2}\right]\leq\mathbb{P}\left[\operatorname{dist}(\hat{Y},\mu)+\operatorname{dist}(\mu,Y^{*})>\alpha_{1}/t+\alpha_{2}/t\right]\leq$
$\displaystyle\mathbb{P}\left[\operatorname{dist}(\hat{Y},\mu)+\sum_{i\in[k]}\frac{w_{i}}{W}\operatorname{dist}(Y_{i},Y^{*})>\alpha_{1}/t+\alpha_{2}/t\right]\leq$
$\displaystyle\mathbb{P}\left[\operatorname{dist}(\hat{Y},\mu)+\alpha_{1}/t>\alpha_{1}/t+\alpha_{2}/t\right]\leq\mathbb{P}\left[\operatorname{dist}(\hat{Y},\mu)>\alpha_{2}/t\right]\leq\beta$

where the first inequality follows from the $t$-approximate triangle
inequality, and the second one follows from the convexity of the
$\operatorname{dist}$ function and the fact that $\mu$ is a convex combination
of $Y_{i}$’s. The third inequality uses the assumption that
$\operatorname{dist}(Y_{i},Y^{*})\leq\alpha_{1}/t$ for all $i\in[k]$. Finally,
the last inequality follows from the assumption that the masking mechanism
$\mathcal{B}$ is $(\alpha_{2}/t,\beta)$ concentrated.

#### Privacy.

Let $D,D^{\prime}$ be two neighbouring datasets and let $\mathcal{M}$ denote
the mechanism specified by Algorithm 1. Note that the statistic $Q$ computed
in Line 4 has sensitivity less than $2/k$. This is because if we modify a
single datapoint, there is at most a single $i$ such that $Y_{i}$ changes.
Hence, the value of $q_{i}$ changes by at most $1$ and for $j\neq i$, the
value of $q_{j}$ changes by less than $1/k$. Thus, $Q$ has sensitivity at most
$2/k$. In particular, since we apply the Truncated Laplace mechanism in Line 7
(see Theorem 2.12), we have

$\mathbb{P}\left[\mathcal{M}(D)=\perp\right]\leq
e^{\varepsilon}\mathbb{P}\left[\mathcal{M}(D^{\prime})=\perp\right]+\delta$
(C.1)

We now show that for any $T\subseteq\mathcal{Y}$, we have

$\displaystyle\mathbb{P}\left[\mathcal{M}(D)\in T\right]$ $\displaystyle\leq
e^{2\varepsilon}\mathbb{P}\left[\mathcal{M}(D^{\prime})\in
T\right]+4e^{\varepsilon}\delta\quad\text{and}$ (C.2)
$\displaystyle\mathbb{P}\left[\mathcal{M}(D)\in T\cup\\{\perp\\}\right]$
$\displaystyle\leq e^{2\varepsilon}\mathbb{P}\left[\mathcal{M}(D^{\prime})\in
T\cup\\{\perp\\}\right]+4e^{\varepsilon}\delta$ (C.3)

which establishes that Algorithm 1 is $(\varepsilon,\delta)$-DP. To this end,
we consider two different cases.

#### Case 1: $Q<0.8$.

In this case,
$\widehat{Q}<0.8+\frac{2}{k\varepsilon}\ln\left(1+\frac{e^{\varepsilon}-1}{2\delta}\right)$
with probability $1$ so $\mathbb{P}\left[\mathcal{M}(D)=\perp\right]=1$. Now,
we verify that Eq. (C.2) and Eq. (C.3) hold. For any $T\subseteq\mathcal{Y}$,
we have $\mathbb{P}\left[\mathcal{M}(D)\in T\right]=0$ so Eq. (C.2) is
trivially satisfied. To check Eq. (C.3) holds, we apply Eq. (C.1) to see that

$\mathbb{P}\left[\mathcal{M}(D)\in
T\cup\\{\perp\\}\right]=\mathbb{P}\left[\mathcal{M}(D)=\perp\right]\leq
e^{\varepsilon}\mathbb{P}\left[\mathcal{M}(D^{\prime})=\perp\right]+\delta\leq
e^{\varepsilon}\mathbb{P}\left[\mathcal{M}(D^{\prime})\in
T\cup\\{\perp\\}\right]+\delta.$

#### Case 2: $Q\geq 0.8$.

Let $Y_{1},\ldots,Y_{k}$ (resp. $Y_{1}^{\prime},\ldots,Y_{k}^{\prime}$) be the
output of the non-private algorithm $\mathcal{A}$ in Line 2 when the dataset
is $D$ (resp. $D^{\prime}$). Further, define let $q_{i},w_{i},\mu$ be as
defined in Algorithm 1 when the input is $D$ and
$q_{i}^{\prime},w_{i}^{\prime},\mu^{\prime}$ be the analogous quantities when
the input is $D^{\prime}$. For notation, let $W=\sum_{i=1}^{k}w_{i}$ and
$W^{\prime}=\sum_{i=1}^{k}W^{\prime}$. Our goal is to show that
$\operatorname{dist}(\mu,\mu^{\prime})\leq 400(r+\phi)/k$ so that the
$(400(r+\phi)/k,\varepsilon,\delta)$-masking mechanism guarantees the output
is $(\varepsilon,\delta)$-indistinguishable.

We proceed by proving a series of claims. The first claim gives a lower bound
on $W$.

######  Claim C.2.

$W\geq k/3$.

###### Proof.

Let $\alpha=|\\{i\in[k]\,:\,q_{i}\geq 0.7\\}|/k$. Then $0.8\leq
Q\leq\alpha+0.7(1-\alpha)$ so $\alpha\geq 1/3$. Finally, $W\geq\alpha k$
because if $q_{i}\geq 0.7$ then $w_{i}=1$. We conclude that $W\geq k/3$. ∎

Let us make one simple observation about the sequence $Y_{1},\ldots,Y_{k}$ and
$Y_{1}^{\prime},\ldots,Y_{k}^{\prime}$. Since only one datapoint has been
changed, it must be the case that $|\\{i\in[k]\,:\,Y_{i}\neq
Y_{i}^{\prime}\\}|\leq 1$. Without loss of generality (perhaps by reindexing),
we assume that $Y_{i}=Y_{i}^{\prime}$ for $2\leq i\leq k$. In addition, we
have the following straightforward claim which follows from the aforementioned
observation.

######  Claim C.3.

For $2\leq i\leq k$, we have $|w_{i}-w_{i}^{\prime}|\leq 10/k$. In addition,
$|W-W^{\prime}|\leq 11$.

###### Proof.

That $|w_{i}-w_{i}^{\prime}|\leq 10/k$ for $2\leq i\leq k$ follows from the
fact that $Y_{i}=Y_{i}^{\prime}$ for $2\leq i\leq k$ (and because Line 8
implies that changing a single datapoint can cause $w_{i}$ can change by at
most $10/k$). Next, $|W-W^{\prime}|\leq|w_{1}-w_{1}^{\prime}|+\sum_{2\leq
i\leq k}|w_{i}-w_{i}^{\prime}|\leq 1+(k-1)\cdot 10/k<11$. ∎

Now, let us proceed to bound $\operatorname{dist}(\mu,\mu^{\prime})$. Using
repeated applications of convexity, we have

$\displaystyle\operatorname{dist}(\mu,\mu^{\prime})$
$\displaystyle=\operatorname{dist}\left(\frac{w_{1}}{W}Y_{1}+\sum_{i=2}^{k}\frac{w_{i}}{W}Y_{i},\frac{w_{1}^{\prime}}{W^{\prime}}Y_{1}^{\prime}+\sum_{j=2}^{k}\frac{w_{j}^{\prime}}{W^{\prime}}Y_{j}\right)$
$\displaystyle\leq\frac{w_{1}}{W}\operatorname{dist}\left(Y_{1},\frac{w_{1}^{\prime}}{W^{\prime}}Y_{1}^{\prime}+\sum_{j=2}^{k}\frac{w_{j}^{\prime}}{W^{\prime}}Y_{j}\right)$
$\displaystyle+\left(1-\frac{w_{1}}{W}\right)\operatorname{dist}\left(\frac{1}{1-w_{1}/W}\sum_{i=2}^{k}\frac{w_{i}}{W}Y_{i},\frac{w_{1}^{\prime}}{W^{\prime}}Y_{1}^{\prime}+\sum_{j=2}^{k}\frac{w_{j}^{\prime}}{W^{\prime}}Y_{j}\right)$
$\displaystyle\leq\frac{w_{1}}{W}\frac{w_{1}^{\prime}}{W^{\prime}}\operatorname{dist}\left(Y_{1},Y_{1}^{\prime}\right)$
(C.4)
$\displaystyle+\frac{w_{1}}{W}\left(1-\frac{w_{1}^{\prime}}{W^{\prime}}\right)\operatorname{dist}\left(Y_{1},\frac{1}{1-w_{1}^{\prime}/W^{\prime}}\sum_{j=2}^{k}\frac{w_{j}^{\prime}}{W^{\prime}}Y_{j}\right)$
(C.5)
$\displaystyle+\left(1-\frac{w_{1}}{W}\right)\frac{w_{1}^{\prime}}{W^{\prime}}\operatorname{dist}\left(\frac{1}{1-w_{1}/W}\sum_{i=2}^{k}\frac{w_{i}}{W}Y_{i},Y_{1}^{\prime}\right)$
(C.6)
$\displaystyle+\left(1-\frac{w_{1}}{W}\right)\left(1-\frac{w_{1}^{\prime}}{W^{\prime}}\right)\operatorname{dist}\left(\frac{1}{1-w_{1}/W}\sum_{i=2}^{k}\frac{w_{i}}{W}Y_{i},\frac{1}{1-w_{1}^{\prime}/W^{\prime}}\sum_{j=2}^{k}\frac{w_{i}^{\prime}}{W^{\prime}}Y_{j}\right),$
(C.7)

where we used convexity in the first argument for the first inequality and
convexity in the second argument for the second inequality. We now bound each
Eq. (C.4), Eq. (C.5), Eq. (C.6), and Eq. (C.7) separately. To do so, we make
use of the following claim.

######  Claim C.4.

Suppose that $w_{i},w_{j}^{\prime}>0$. Then
$\operatorname{dist}(Y_{i},Y_{j}^{\prime})\leq 2r$.

###### Proof.

Let $S_{i}=\\{s\in[k]\,:\,\operatorname{dist}(Y_{i},Y_{s})\leq r/t\\}$ and
$S_{j}^{\prime}=\\{s\in[k]\,:\,\operatorname{dist}(Y_{j}^{\prime},Y_{s}^{\prime})\leq
r/t\\}$. Since $w_{i},w_{j}^{\prime}>0$, we must have $|S_{i}|\geq 0.6k$ and
$|S_{j}^{\prime}|\geq 0.6k$. Hence, $|S_{i}\cap
S_{j}^{\prime}|=|S_{i}|+|S_{j}^{\prime}|-|S_{i}\cup S_{j}^{\prime}|\geq
0.2k\geq 2$ since $k\geq 10$. So there exists $s\in S_{i}\cap S_{j}^{\prime}$
such that $2\leq s\leq k$ which implies that $Y_{s}=Y_{s}^{\prime}$. In
particular, for this $s$, we have $\operatorname{dist}(Y_{i},Y_{s})\leq r/t$
and $\operatorname{dist}(Y_{j}^{\prime},Y_{s})\leq r/t$. Recalling that
$\operatorname{dist}$ satisfies a $t$-approximate $r$-restricted triangle
inequality, we have $\operatorname{dist}(Y_{i},Y_{j}^{\prime})\leq
t(\operatorname{dist}(Y_{i},Y_{s})+\operatorname{dist}(Y_{j}^{\prime},Y_{s}))\leq
2r$. ∎

By Claim C.4, we have that Eq. (C.4) is bounded by

$\frac{w_{1}w_{1}^{\prime}}{WW^{\prime}}\operatorname{dist}(Y_{1},Y_{1}^{\prime})\leq\frac{w_{1}w_{1}^{\prime}}{WW^{\prime}}2r$
(C.8)

We now bound Eq. (C.5) which is very similar to the bound of Eq. (C.4) we just
obtained. Using convexity, we have

$\displaystyle\frac{w_{1}}{W}\left(1-\frac{w_{1}^{\prime}}{W^{\prime}}\right)\operatorname{dist}\left(Y_{1},\frac{1}{1-w_{1}^{\prime}/W^{\prime}}\sum_{j=2}^{k}\frac{w_{j}^{\prime}}{W^{\prime}}Y_{j}\right)$
$\displaystyle\leq\frac{w_{1}}{W}\sum_{j=2}^{k}\frac{w_{j}^{\prime}}{W^{\prime}}\operatorname{dist}(Y_{1},Y_{j})$
(C.9)
$\displaystyle\leq\frac{w_{1}}{W}\sum_{j=2}^{k}\frac{w_{j}^{\prime}}{W^{\prime}}\cdot
2r,$

where we applied Claim C.4 for the last inequality. Similarly, we have

$\displaystyle\left(1-\frac{w_{1}}{W}\right)\frac{w_{1}^{\prime}}{W^{\prime}}\operatorname{dist}\left(\frac{1}{1-w_{1}/W}\sum_{i=2}^{k}\frac{w_{i}}{W}Y_{i},Y_{1}^{\prime}\right)\leq\frac{w_{1}^{\prime}}{W^{\prime}}\sum_{i=2}^{k}\frac{w_{j}}{W}2r.$
(C.10)

Note that Eq. (C.8), Eq. (C.9), Eq. (C.10) shows that the sum of Eq. (C.4),
Eq. (C.5), Eq. (C.6) is bounded above by $\frac{2r}{W}+\frac{2r}{W^{\prime}}$.

Finally, we bound Eq. (C.7). We require the following claim whose proof is
similar to that of Claim C.4.

######  Claim C.5.

Let $S=\\{2\leq i\leq k\,:\,\max(w_{i},w_{i}^{\prime})>0\\}$. Then for all
$i,j\in S$, we have $\operatorname{dist}(Y_{i},Y_{j})\leq 2r$.

###### Proof.

Let $S_{i}=\\{s\in[k]\,:\,\operatorname{dist}(Y_{i},Y_{s})\leq r/t\\}$ and
$S_{j}=\\{s\in[k]\,:\,\operatorname{dist}(Y_{j},Y_{s})\leq r/t\\}$. Since
$\max(w_{i},w_{i}^{\prime}),\max(w_{j},w_{j}^{\prime})>0$, we must have
$|S_{i}|\geq 0.6k-1$ and $|S_{j}|\geq 0.6k-1$. Hence, $|S_{i}\cap
S_{j}|=|S_{i}|+|S_{j}|-|S_{i}\cup S_{j}|\geq 0.2k-2>1$ so $|S_{i}\cap
S_{j}|\geq 2$. So there exists $s\in S_{i}\cap S_{j}$ such that $s\geq 2$. In
particular, for this $s$, we have $\operatorname{dist}(Y_{i},Y_{s})\leq r/t$
and $\operatorname{dist}(Y_{j},Y_{s})\leq r/t$. Recalling that
$\operatorname{dist}$ satisfies a $t$-approximate $r$-restricted triangle
inequality, we have $\operatorname{dist}(Y_{i},Y_{j})\leq
t(\operatorname{dist}(Y_{i},Y_{s})+\operatorname{dist}(Y_{j},Y_{s}))\leq 2r$.
∎

Note that $\sum_{i=2}^{k}\frac{w_{i}}{W}Y_{i}=\sum_{i\in
S}\frac{w_{i}}{W}Y_{i}$ and
$\sum_{i=2}^{k}\frac{w_{i}^{\prime}}{W^{\prime}}Y_{i}=\sum_{i\in
S}\frac{w_{i}^{\prime}}{W^{\prime}}Y_{i}$. Now, we have

$\displaystyle\operatorname{dist}$
$\displaystyle\left(\frac{1}{1-w_{1}/W}\sum_{i\in
S}\frac{w_{i}}{W}Y_{i},\frac{1}{1-w_{1}^{\prime}/W^{\prime}}\sum_{i\in
S}\frac{w_{i}^{\prime}}{W^{\prime}}Y_{i}\right)$
$\displaystyle\leq\left({\phi}+\max_{i,j\in
S}\operatorname{dist}(Y_{i},Y_{j})\right)\sum_{i\in
S}\left|\frac{w_{i}/W}{1-w_{1}/W}-\frac{w_{i}^{\prime}/W^{\prime}}{1-w_{1}^{\prime}/W^{\prime}}\right|.$

Where we used the fact that $\operatorname{dist}$ satisfies locality in
Definition 3.1. Note that $\operatorname{dist}(Y_{i},Y_{j})\leq 2r$ for all
$i,j\in S$. Also, by the Claim C.6 we see that the right hand side is bounded
by $100/k$. Therefore, the term in Eq. (C.7) is upper bounded by
$100(2r+\phi)/k$. We now state and prove Claim C.6.

######  Claim C.6.

$\sum_{i\in
S}\left|\frac{w_{i}/W}{1-w_{1}/W}-\frac{w_{i}^{\prime}/W^{\prime}}{1-w_{1}^{\prime}/W^{\prime}}\right|\leq
100/k$.

###### Proof.

$\displaystyle\sum_{i\in
S}\left|\frac{w_{i}/W}{1-w_{1}/W}-\frac{w_{i}^{\prime}/W^{\prime}}{1-w_{1}^{\prime}/W^{\prime}}\right|=\sum_{i\in
S}\left|\frac{w_{i}}{W-w_{1}}-\frac{w_{i}^{\prime}}{W^{\prime}-w_{1}^{\prime}}\right|$
$\displaystyle\leq\sum_{i\in
S}\left(\left|\frac{w_{i}}{W-w_{1}}-\frac{w_{i}}{W^{\prime}-w_{1}^{\prime}}\right|+\left|\frac{w_{i}}{W^{\prime}-w_{1}^{\prime}}-\frac{w_{i}^{\prime}}{W^{\prime}-w_{1}^{\prime}}\right|\right)$
$\displaystyle\leq\sum_{i\in
S}\left(\frac{w_{i}}{W}\left|\frac{(W^{\prime}-W)+(w_{1}-w_{1}^{\prime})}{(1-w_{1}/W)(W^{\prime}-w_{1}^{\prime})}\right|+|w_{i}-w_{i}^{\prime}|\left|\frac{1}{W^{\prime}-w_{1}^{\prime}}\right|\right)$
$\displaystyle\leq\sum_{i\in
S}\left(\frac{w_{i}}{W}\left|\frac{11+2}{(1-3/k)(k/3-12)}\right|+|10/k|\frac{1}{k/3-12}\right)$
$\displaystyle\leq\left(\frac{39}{0.9(k-36)}\right)+\frac{30}{k-36}\leq\frac{74}{k-36}\leq
100/k.$

In the second last line, we used that $|1-w_{1}/W|=1-w_{1}/W\geq 1-1/W\geq
1-3/k$ by Claim C.2 (provided $k>3$) and
$|W^{\prime}-w_{1}^{\prime}|=W^{\prime}-w_{1}^{\prime}\geq
W-|W-W^{\prime}|-w_{1}^{\prime}\geq k/3-12$ by Claim C.2 and Claim C.3
(provided $k>36$). For the last line, we used that $\sum_{i\in S}w_{i}/W\leq
1$, $\sum{i\in S}1/k\leq 1$, and $1-3/k\geq 0.9$ for $k\geq 30$. Finally, the
last inequality is true provided $k\geq 140$.

We have now bounded Eq. (C.8), Eq. (C.9), Eq. (C.10), and Eq. (C.7). The sum
of these terms is at most $\frac{2r}{W}+\frac{2r}{W^{\prime}}+100(2r+\phi)/k$.
Using the Claim C.2 again we can upper bound this sum by
$(212r+100\phi)/k<400(r+\phi)/k$. ∎

We now verify that the algorithm is
$(2\varepsilon,{4e^{\varepsilon}\delta})$-DP in the present case. Since
$\operatorname{dist}(\mu,\mu^{\prime})\leq 400(r+\phi)/k$ and $\mathcal{B}$ is
a $(400(r+\phi)/k,\varepsilon,\delta)$-masking mechanism, we know that
$\mathcal{B}(\mu),\mathcal{B}(\mu^{\prime})$ are
$(\varepsilon,\delta)$-indistinguishable. Let $T\subseteq\mathcal{Y}$. Then

$\displaystyle\mathbb{P}\left[\mathcal{M}(D)\in T\right]$
$\displaystyle=\mathbb{P}\left[\mathcal{M}(D)\neq\perp\right]\mathbb{P}\left[\mathcal{B}(\mu)\in
T\right]$
$\displaystyle\leq(e^{\varepsilon}\mathbb{P}\left[\mathcal{M}(D^{\prime})\neq\perp\right]+\delta)(e^{\varepsilon}\mathbb{P}\left[\mathcal{B}(\mu^{\prime})\in
T\right]+\delta)$
$\displaystyle=e^{2\varepsilon}\mathbb{P}\left[\mathcal{M}(D^{\prime})\neq\perp\right]\mathbb{P}\left[\mathcal{B}(\mu^{\prime})\in
T\right]$
$\displaystyle+e^{\varepsilon}\delta(\mathbb{P}\left[\mathcal{M}(D^{\prime})\neq\perp\right]+\mathbb{P}\left[\mathcal{B}(\mu^{\prime})\in
T\right])+\delta^{2}$ $\displaystyle\leq
e^{2\varepsilon}\mathbb{P}\left[\mathcal{M}(D^{\prime})\in
T\right]+3e^{\varepsilon}\delta,$

which establishes Eq. (C.2). Finally,

$\displaystyle\mathbb{P}\left[\mathcal{M}(D)\in T\cup\\{\perp\\}\right]$
$\displaystyle=\mathbb{P}\left[\mathcal{M}(D)=\perp\right]+\mathbb{P}\left[\mathcal{M}(D)\in
T\right]$
$\displaystyle\leq(e^{\varepsilon}\mathbb{P}\left[\mathcal{M}(D^{\prime})=\perp\right]+\delta)+(e^{2\varepsilon}\mathbb{P}\left[\mathcal{M}(D^{\prime})\in
T\right]+3e^{\varepsilon}\delta)$ $\displaystyle\leq
e^{2\varepsilon}\mathbb{P}\left[\mathcal{M}(D^{\prime})\in
T\cup\\{\perp\\}\right]+4e^{\varepsilon}\delta,$

which establishes Eq. (C.3). ∎

##  Appendix D Missing proofs from Section 5

###  D.1 Proof of 5.1

5.1 follows easily from the following lemma.

######  Lemma D.1.

Suppose that

$\max\\{\|\Sigma_{1}^{-1/2}\Sigma_{2}\Sigma_{1}^{-1/2}-I_{d}\|,\|\Sigma_{2}^{-1/2}\Sigma_{1}\Sigma_{2}^{-1/2}-I_{d}\|\\}\leq\gamma\leq
1/2.$

Then for every $\delta_{0}>0$, $\Sigma_{1}^{1/2}(I+\eta G)$ and
$\Sigma_{2}^{1/2}(I+\eta G)$ are
$(\varepsilon_{0}(\delta_{0}),\delta_{0})$-indistinguishable where

$\varepsilon_{0}(\delta_{0})=\frac{\gamma^{2}d}{2}(d+1/\eta^{2})+2d\gamma\sqrt{\ln(2/\delta_{0})}+2\gamma\ln(2/\delta_{0})+3\gamma\sqrt{d}\sqrt{\ln(2/\delta_{0})}/\eta.$
(D.1)

Indeed, given D.1, 5.1 follows by setting $\gamma$ appropriately so that each
term in Eq. (D.1) is at most $\varepsilon/4$ and post-processing, i.e. since
$\Sigma_{1}^{1/2}(I+\eta G)$ and $\Sigma_{2}^{1/2}(I+\eta G)$ are
$(\varepsilon,\delta)$-indistinguishable then so are $\Sigma_{1}^{1/2}(I+\eta
G)(I+\eta G)^{\top}\Sigma_{1}^{1/2}$ and $\Sigma_{2}^{1/2}(I+\eta G)(I+\eta
G)^{\top}\Sigma_{2}^{1/2}$.

First, we require the following technical lemma which asserts that for
Gaussians
$\mathcal{N}(\mu,\Sigma),\mathcal{N}(\widehat{\mu},\widehat{\Sigma})$, if
$g\sim\mathcal{N}(0,I_{d})$ then the random variable
$\mathcal{L}_{{\mathcal{N}(\mu,\Sigma)}\parallel{\mathcal{N}(\widehat{\mu},\widehat{\Sigma})}}(\Sigma^{1/2}g+\mu)$
concentrates tightly around
$D_{\mathrm{KL}}\left({\mathcal{N}(\mu,\Sigma)}\parallel{\mathcal{N}(\widehat{\mu},\widehat{\Sigma})}\right)$.
The plan is to apply the lemma with $\mu\in\mathbb{R}^{d^{2}}$ (resp.
$\widehat{\mu}$) is $\Sigma_{1}^{1/2}$ (resp. $\Sigma_{2}^{1/2}$) but viewed
as a vector of length $d^{2}$ and $\Sigma$ (resp. $\widehat{\Sigma}$) is the
covariance matrix of (the vector) $\Sigma_{1}^{1/2}(I+\eta G)$ (resp.
$\Sigma_{2}^{1/2}(I+\eta G)$).

######  Lemma D.2.

Let $\mu,\widehat{\mu}\in\mathbb{R}^{d}$ and
$\Sigma,\widehat{\Sigma}\in\mathbb{S}_{>0}^{d}$. Let
$g\sim\mathcal{N}(0,I_{d})$. Then

$\displaystyle\mathcal{L}_{{\mathcal{N}(\mu,\Sigma)}\parallel{\mathcal{N}(\widehat{\mu},\widehat{\Sigma})}}(\Sigma^{1/2}g+\mu)$
$\displaystyle\leq
D_{\mathrm{KL}}\left({\mathcal{N}(\mu,\Sigma)}\parallel{\mathcal{N}(\widehat{\mu},\widehat{\Sigma})}\right)$
$\displaystyle+2\|\Sigma^{1/2}\widehat{\Sigma}^{-1}\Sigma^{1/2}-I\|_{F}\sqrt{\ln(2/\delta)}+2\|\Sigma^{1/2}\widehat{\Sigma}^{-1}\Sigma^{1/2}-I\|\ln(2/\delta)$
$\displaystyle+\|\Sigma^{1/2}\widehat{\Sigma}^{-1}\Sigma^{1/2}\|\cdot\|\Sigma^{-1/2}(\widehat{\mu}-\mu)\|_{2}\sqrt{2\ln(2/\delta)}$

with probability at least $1-\delta$.

The proof of D.2 can be found in Subsection D.2.

We use D.2 to show that if
$\|\Sigma_{1}^{-1/2}\Sigma_{2}\Sigma_{1}^{-1/2}-I\|$ and
$\|\Sigma_{2}^{-1/2}\Sigma_{1}\Sigma_{2}^{-1/2}-I\|$ are small then
$\Sigma_{1}^{1/2}(I_{d}+\eta G)$ and $\Sigma_{2}^{1/2}(I_{d}+\eta G)$ are
indistinguishable. As a consequence, the matrices $\Sigma_{1}^{1/2}(I_{d}+\eta
G)(I_{d}+\eta G)^{\top}\Sigma_{1}^{1/2}$ and $\Sigma_{2}^{1/2}(I_{d}+\eta
G)(I_{d}+\eta G)^{\top}\Sigma_{2}^{1/2}$ are also indistinguishable.

Before we prove D.1, it is convenient to set up some notation. For a matrix
$B=\begin{bmatrix}b_{1}\cdots b_{d}\end{bmatrix}\in\mathbb{R}^{d\times d}$, we
write $\operatorname{vec}(B)$ to denote a vector of length $d^{2}$ with the
columns of $B$ stacked together. More precisely,

$\operatorname{vec}(B)=\begin{bmatrix}b_{1}\\\ \vdots\\\
b_{d}\end{bmatrix}\in\mathbb{R}^{d^{2}}.$

Next, for a matrix $A=\begin{bmatrix}a_{1}\cdots
a_{d}\end{bmatrix}\in\mathbb{R}^{k\times d}$, we let $\operatorname{bd}(A)$ be
the block-diagonal matrix with $d$ copies of $A$ along its diagonal, i.e.

$\operatorname{bd}(A)=\begin{bmatrix}A&0&\ldots&0\\\ 0&A&\ldots&0\\\
\vdots&\vdots&\ddots&\vdots\\\
0&0&\ldots&A\end{bmatrix}\in\mathbb{R}^{kd\times d^{2}}.$

Note that for a size-compatible matrix
$\operatorname{vec}(AB)=\operatorname{bd}(A)\operatorname{vec}(B)$. Indeed,
$AB=\begin{bmatrix}Ab_{1}\cdots Ab_{d}\end{bmatrix}$ so

$\underbrace{\begin{bmatrix}Ab_{1}\\\ \vdots\\\
Ab_{d}\end{bmatrix}}_{\operatorname{vec}(AB)}=\underbrace{\begin{bmatrix}A&0&\ldots&0\\\
0&A&\ldots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\ldots&A\end{bmatrix}}_{\operatorname{bd}(A)}\underbrace{\begin{bmatrix}b_{1}\\\
\vdots\\\ b_{d}\end{bmatrix}}_{\operatorname{vec}(B)}.$

We require a few calculations before we can prove D.1.

######  Claim D.3.

Suppose that
$\|\Sigma_{1}^{1/2}\Sigma_{2}^{-1}\Sigma_{1}^{1/2}-I_{d}\|\leq\gamma$. Then

$\|\operatorname{bd}(\Sigma_{1})^{1/2}\operatorname{bd}(\Sigma_{2})^{-1}\operatorname{bd}(\Sigma_{1})^{1/2}-I_{d^{2}}\|_{F}^{2}=d\cdot\|\Sigma_{1}^{-1/2}\Sigma_{2}\Sigma_{1}^{-1/2}-I_{d}\|_{F}^{2}\leq
d^{2}\gamma^{2}.$

###### Proof.

The equality is because the matrices
$\operatorname{bd}(\Sigma_{1}),\operatorname{bd}(\Sigma_{2}),I_{d^{2}}$ are
all block diagonal with copies of $\Sigma_{1},\Sigma_{2},I_{d}$ along its
diagonal.

For the inequality, we have
$\|\Sigma_{1}^{1/2}\Sigma_{2}^{-1}\Sigma_{1}^{1/2}-I_{d}\|\leq\gamma$ which
implies
$\|\Sigma_{1}^{1/2}\Sigma_{2}^{-1}\Sigma_{1}^{1/2}-I_{d}\|_{F}^{2}\leq\gamma^{2}\cdot
d$. Multiplying both sides by $d$ gives the result. ∎

######  Claim D.4.

Suppose that
$\|\Sigma_{1}^{1/2}\Sigma_{2}^{-1}\Sigma_{1}^{1/2}-I_{d}\|\leq\gamma$. Then
$\|\eta^{-1}\operatorname{bd}(\Sigma_{2})^{-1/2}(\operatorname{vec}(\Sigma_{2}^{1/2})-\operatorname{vec}(\Sigma_{1}^{1/2}))\|_{2}^{2}\leq\gamma^{2}d/\eta^{2}$.

###### Proof.

Note that
$\|\operatorname{bd}(\Sigma_{2})^{-1/2}(\operatorname{vec}(\Sigma_{2}^{1/2})-\operatorname{vec}(\Sigma_{1}^{1/2}))\|_{2}^{2}=\|I-\Sigma_{2}^{-1/2}\Sigma_{1}^{1/2}\|_{F}^{2}$.
Let $\lambda_{1},\ldots,\lambda_{d}$ be the (positive) eigenvalues of
$\Sigma_{1}^{1/2}\Sigma_{2}^{-1}\Sigma_{1}^{1/2}$.555 The eigenvalues of
$\Sigma_{2}^{-1/2}\Sigma_{1}\Sigma_{2}^{-1/2}-I_{d}$ lie in the interval
$[-1/2,1/2]$ so the eigenvalues of
$\Sigma_{2}^{-1/2}\Sigma_{1}\Sigma_{2}^{-1/2}$ lie in the interval
$[1/2,3/2]$. In particular, the eigenvalues are positive. Then
$\sqrt{\lambda_{1}},\ldots,\sqrt{\lambda_{d}}$ are the singular values of
$\Sigma_{2}^{-1/2}\Sigma_{1}^{1/2}$ and
$1-\sqrt{\lambda_{1}},\ldots,1-\sqrt{\lambda_{d}}$ are the singular values of
$I_{d}-\Sigma_{2}^{-1/2}\Sigma_{1}^{1/2}$. Thus,

$\displaystyle\|I_{d}-\Sigma_{2}^{-1/2}\Sigma_{1}^{1/2}\|_{F}^{2}$
$\displaystyle=\sum_{i=1}^{d}(1-\sqrt{\lambda_{i}})^{2}$
$\displaystyle=\sum_{i=1}^{d}1-2\sqrt{\lambda_{i}}+\lambda_{i}$
$\displaystyle\leq\sum_{i=1}^{d}1-2\sqrt{1-\gamma}+(1+\gamma)$
$\displaystyle\leq\sum_{i=1}^{d}1-2(1-\gamma/2-\gamma^{2}/2)+(1+\gamma)$
$\displaystyle=\sum_{i=1}^{d}\gamma^{2}=\gamma^{2}d,$

where in the third line we used that $\lambda_{i}\in[1-\gamma,1+\gamma]$ for
all $i\in[d]$ and in the fourth line, we used the assumption that $\gamma\leq
1/2$ and Claim B.2 to bound $\sqrt{1-\gamma}\leq 1-\gamma/2-\gamma^{2}/2$. The
lemma follows by dividing both sides by $\eta^{2}$. ∎

######  Claim D.5.

Suppose that
$\|\Sigma_{1}^{1/2}\Sigma_{2}^{-1}\Sigma_{1}^{1/2}-I_{d}\|\leq\gamma\leq 1/2$.
Then

$D_{\mathrm{KL}}\left({\mathcal{N}(\operatorname{vec}(\Sigma_{1}^{1/2}),\eta^{2}\operatorname{bd}(\Sigma_{1}))}\parallel{\mathcal{N}(\operatorname{vec}(\Sigma_{2}^{1/2}),\eta^{2}\operatorname{bd}(\Sigma_{2}))}\right)\leq\frac{\gamma^{2}d}{2}(d+1/\eta^{2}).$

###### Proof.

Follows by plugging in the bounds from Claim D.3 and Claim D.4 into 2.1. ∎

###### Proof of D.1.

It is convenient to view $\Sigma_{i}^{1/2}(I_{d}+\eta G)$ as
$d^{2}$-dimensional vectors. Indeed, we have
$\operatorname{vec}(\Sigma_{i}^{1/2}(I_{d}+\eta
G))=\operatorname{vec}(\Sigma_{i}^{1/2})+\eta\operatorname{bd}(\Sigma_{i}^{1/2})\operatorname{vec}(G)$
which is a random $d^{2}$-dimensional vector distributed as
$\mathcal{N}(\operatorname{vec}(\Sigma_{i}^{1/2}),\eta^{2}\operatorname{bd}(\Sigma_{i}))$.
For notation, let us write
$\mathcal{D}_{i}=\mathcal{N}(\operatorname{vec}(\Sigma_{i}^{1/2}),\eta^{2}\operatorname{bd}(\Sigma_{i}))$.
We will use D.2 to bound

$\mathcal{L}\coloneqq\mathcal{L}_{{\mathcal{D}_{1}}\parallel{\mathcal{D}_{2}}}(\eta\operatorname{bd}(\Sigma_{1})^{1/2}\operatorname{vec}(G)+\operatorname{vec}(\Sigma_{1})^{1/2}).$

An identical argument holds with $\Sigma_{1},\Sigma_{2}$ reversed and so we
may apply 2.10.

By D.2, we have that

$\displaystyle\mathcal{L}$ $\displaystyle\leq
D_{\mathrm{KL}}\left({\mathcal{D}_{1}}\parallel{\mathcal{D}_{2}}\right)$ (D.2)
$\displaystyle+2\|\operatorname{bd}(\Sigma_{1})^{1/2}\operatorname{bd}(\Sigma_{2})^{-1}\operatorname{bd}(\Sigma_{1})^{1/2}-I\|_{F}\sqrt{\ln(2/\delta)}$
$\displaystyle+2\|\operatorname{bd}(\Sigma_{1})^{1/2}\operatorname{bd}(\Sigma_{2})^{-1}\operatorname{bd}(\Sigma_{1})^{1/2}-I_{d^{2}}\|\ln(2/\delta)$
$\displaystyle+\|\operatorname{bd}(\Sigma_{1})^{1/2}\operatorname{bd}(\Sigma_{2})^{-1}\operatorname{bd}(\Sigma_{1})^{1/2}\|\cdot\|\eta^{-1}\operatorname{bd}(\Sigma_{1})^{-1/2}(\operatorname{vec}(\Sigma_{2})-\operatorname{vec}(\Sigma_{1}))\|_{2}\sqrt{2\ln(2/\delta)},$

with probability at least $1-\delta$. We now bound each term in Eq. (D.2).

The first term in Eq. (D.2) is bounded by Claim D.5 and we have

$D_{\mathrm{KL}}\left({\mathcal{D}_{1}}\parallel{\mathcal{D}_{2}}\right)\leq\frac{\gamma^{2}d}{2}(d+1/\eta^{2}).$
(D.3)

The second term is bounded by Claim D.3 and we have

$\|\operatorname{bd}(\Sigma_{1})^{1/2}\operatorname{bd}(\Sigma_{2})^{-1}\operatorname{bd}(\Sigma_{1})^{1/2}-I_{d^{2}}\|_{F}\leq\gamma
d.$ (D.4)

For the third term, note that
$\operatorname{bd}(\Sigma_{1})^{1/2}\operatorname{bd}(\Sigma_{2})^{-1}\operatorname{bd}(\Sigma_{1})^{1/2}$
and $\Sigma_{1}^{1/2}\Sigma_{2}^{-1}\Sigma_{1}$ have the same set of
eigenvalues but with different multiplicities (by a factor of $d$). Hence,

$\|\operatorname{bd}(\Sigma_{1})^{1/2}\operatorname{bd}(\Sigma_{2})^{-1}\operatorname{bd}(\Sigma_{1})^{1/2}-I_{d^{2}}\|\leq\gamma.$
(D.5)

Finally, for the last term, we use that
$\|\operatorname{bd}(\Sigma_{1})^{1/2}\operatorname{bd}(\Sigma_{2})^{-1}\operatorname{bd}(\Sigma_{1})^{1/2}\|\leq
1+\gamma\leq 3/2$ and Claim D.4 to obtain

$\|\operatorname{bd}(\Sigma_{1})^{1/2}\operatorname{bd}(\Sigma_{2})^{-1}\operatorname{bd}(\Sigma_{1})^{1/2}\|\cdot\|\eta^{-1}\operatorname{bd}(\Sigma_{1})^{-1/2}(\operatorname{vec}(\Sigma_{2})-\operatorname{vec}(\Sigma_{1}))\|_{2}<3\gamma\sqrt{d}/2\eta.$
(D.6)

Plugging Eq. (D.3) through Eq. (D.6) into Eq. (D.2) gives

$\mathcal{L}\leq\frac{\gamma^{2}d}{2}(d+1/\eta^{2})+2d\gamma\sqrt{\ln(2/\delta)}+2\gamma\ln(2/\delta)+3\gamma\sqrt{d}\sqrt{\ln(2/\delta)}/\eta,$

with probability at least $1-\delta$. ∎

###  D.2 Proof of D.2

###### Proof of D.2.

Let
$\mathcal{L}=\mathcal{L}_{{\mathcal{N}(\mu,\Sigma)}\parallel{\mathcal{N}(\widehat{\mu},\widehat{\Sigma})}}(\Sigma^{1/2}g+\mu)$.
From Claim 2.11, we have that

$\mathcal{L}=-\frac{1}{2}\ln(\det(\widehat{\Sigma}^{-1}\Sigma))-\frac{1}{2}g^{\top}g+\frac{1}{2}(g-\Sigma^{-1/2}(\widehat{\mu}-\mu))^{\top}\Sigma^{1/2}\widehat{\Sigma}^{-1}\Sigma^{1/2}(g-\Sigma^{-1/2}(\widehat{\mu}-\mu)).$
(D.7)

Expanding the last term, we have

$\mathcal{L}=-\frac{1}{2}\ln(\det(\widehat{\Sigma}^{-1}\Sigma))+\frac{1}{2}g^{\top}(\Sigma^{1/2}\widehat{\Sigma}^{-1}\Sigma^{1/2}-I)g-(\widehat{\mu}-\mu)^{\top}\widehat{\Sigma}^{-1}\Sigma^{1/2}g+\frac{1}{2}(\widehat{\mu}-\mu)^{\top}\widehat{\Sigma}^{-1}(\widehat{\mu}-\mu).$
(D.8)

By 2.4, we have that

$\displaystyle g^{\top}(\Sigma^{1/2}\widehat{\Sigma}^{-1}\Sigma^{1/2}-I)g$
$\displaystyle\leq\operatorname{tr}(\Sigma^{1/2}\widehat{\Sigma}^{-1}\Sigma^{1/2}-I)+4\|\Sigma^{1/2}\widehat{\Sigma}^{-1}\Sigma^{1/2}-I\|_{F}\sqrt{\ln(2/\delta)}$
(D.9)
$\displaystyle+4\|\Sigma^{1/2}\widehat{\Sigma}^{-1}\Sigma^{1/2}-I\|\ln(2/\delta),$

with probability $1-\delta/2$. Next,
$(\widehat{\mu}-\mu)^{\top}\widehat{\Sigma}^{-1}\Sigma^{1/2}g$ is a mean-zero
Gaussian random variable with variance
$\|\Sigma^{1/2}\widehat{\Sigma}^{-1}(\widehat{\mu}-\mu)\|_{2}^{2}\leq\|\Sigma^{1/2}\widehat{\Sigma}^{-1}\Sigma^{1/2}\|^{2}\|\Sigma^{-1/2}(\widehat{\mu}-\mu)\|_{2}^{2}$.
Hence, using the tails of a Gaussian random variable (2.3), we have

$(\widehat{\mu}-\mu)^{\top}\widehat{\Sigma}^{-1}\Sigma^{1/2}g\leq\|\Sigma^{1/2}\widehat{\Sigma}^{-1}\Sigma^{1/2}\|\|\Sigma^{-1/2}(\widehat{\mu}-\mu)\|_{2}\sqrt{2\ln(2/\delta)}$
(D.10)

with probability $1-\delta/2$. Plugging Eq. (D.9) and Eq. (D.10) along with
the formula for the KL divergence between two multivariate Normal
distributions (2.1), we have

$\displaystyle\mathcal{L}$ $\displaystyle\leq
D_{\mathrm{KL}}\left({\mathcal{N}(\mu,\Sigma)}\parallel{\mathcal{N}(\widehat{\mu},\widehat{\Sigma})}\right)+2\|\Sigma^{1/2}\widehat{\Sigma}^{-1}\Sigma^{1/2}-I\|_{F}\sqrt{\ln(2/\delta)}$
(D.11)
$\displaystyle+2\|\Sigma^{1/2}\widehat{\Sigma}^{-1}\Sigma^{1/2}-I\|\ln(2/\delta)+\|\Sigma^{1/2}\widehat{\Sigma}^{-1}\Sigma^{1/2}\|\|\Sigma^{-1/2}(\widehat{\mu}-\mu)\|_{2}\sqrt{2\ln(2/\delta)},$

with probability $1-\delta$. ∎