CoinPress: Practical Private Mean and Covariance Estimation

Introduction

One of the most basic problems in statistics and machine learning is to estimate the mean and covariance of a distribution based on i.i.d. samples. Not only are these some of the most basic summary statistics one could want for a real-valued distribution, but they are also building blocks for more sophisticated statistical estimation tasks like linear regression and stochastic convex optimization.

The optimal solutions to these problems are folklore—simply output the empirical mean and covariance of the samples. However, this solution is not suitable when the samples consist of sensitive, private information belonging to individuals, as it has been shown repeatedly that even releasing just the empirical mean can reveal this sensitive information [DN03, HSR+08, BUV14, DSS+15, DSSU17]. Thus, we need estimators that are not only accurate with respect to the underlying distribution, but also protect the privacy of the individuals represented in the sample.

The most widely accepted solution to individual privacy in statistics and machine learning is differential privacy (DP) [DMNS06], which provides a strong guarantee of individual privacy by ensuring that no individual has a significant influence on the learned parameters. A large body of work now shows that, in principle, nearly every statistical task can be solved privately, and differential privacy is now being deployed by Apple [Dif17], Google [EPK14, BEM+17], Microsoft [DKY17], and the US Census Bureau [DLS+17].

Differential privacy requires adding random noise to some stage of the estimation procedure, and this noise might increase the error of the final estimate. Typically, the amount of noise vanishes as the sample size $n$ grows, and one can often show that as $n\to\infty$, the additional error due to privacy vanishes faster than the sampling error of the estimator, making differential privacy highly practical for large samples.

However, differential privacy is often difficult to achieve for small datasets, or when the dataset is large, but we want to restrict attention to some small subpopulation within the data. Thus, a recent trend has been to focus on simple, widely used estimation tasks, and design estimators with good concrete performance at small samples sizes. Most relevant to our work, Karwa and Vadhan [KV18] and Du, Foot, Moniot, Bray, and Groce [DFM+20] give practical mean and variance estimators for univariate Gaussian data. However, as we show, these methods do not scale well to the more challenging multivariate setting.

In this work we give simple, practical estimators for the mean and covariance of multivariate sub-Gaussian data. We call our method CoinPress, for COnfidence-INterval-based PRivate EStimation Strategy. We validate our estimators theoretically and empirically. On the theoretical side, we show that our estimators match the state-of-the-art asymptotic bounds for sub-Gaussian mean and covariance estimation [KLSU19]. On the empirical side, we give an extensive evaluation with synthetic data, as well as a demonstration on a real-world dataset. We show that our estimators have error comparable to that of the non-private empirical mean and covariance at small sample sizes. See Figure 1 for one representative example of our algorithm’s performance. Our mean estimator also improves over the state-of-the-art method of Du et al. [DFM+20], which was developed for univariate data but can be applied coordinate-wise to estimate multivariate data. We highlight a few other important features of our methods:

First, like many differentially private estimators, our method requires the user to input some a priori knowledge of the data. For mean estimation, we require the mean lives in a specified ball of radius $R$, and for covariance estimation we require that the covariance matrix can be sandwiched spectrally between $A$ and $KA$ for some matrix $A$. Some a priori boundedness is necessary for algorithms like ours that satisfy concentrated DP [DR16, BS16, Mir17, BDRS18], or satisfy pure DP. Under pure or concentrated DP, the dependence on $R$ and $K$ must be polylogarithmic [KV18, BKSW19]. One can allow $R=\infty$ for mean estimation under $(\varepsilon,\delta)$-DP with $\delta>0$ [KV18], although the resulting algorithm has poor concrete performance even for univariate data. It is an open question whether one can allow $K=\infty$ for covariance estimation even under $(\varepsilon,\delta)$-DP. We show that our estimator is practical when these parameters are taken to be extremely large, meaning the user only needs a very weak prior.

Second, for simplicity, we describe and evaluate our methods primarily with Gaussian data. However, the only feature of Gaussian data that is relevant for our methods is a strong bound on the tails of the distribution, and, by definition, these bounds hold for any sub-Gaussian distribution. Moreover, using experiments with both heavier-tailed synthetic data and with real-world data, we demonstrate that our method remains useful even when the data is not truly Gaussian. Note that some restriction on the details of the data is necessary, at least in the worst-case, as [KSU20] showed that the minimax optimal error is highly sensitive to the rate of decay of the distribution’s tails.

Approach. At a high-level, our estimators work by iteratively refining an estimate for the parameters, inspired by [KLSU19]. For mean estimation, we start with some (potentially very large) ball $B_{1}$ of radius $R_{1}$ that we know contains most of the mass of the probability distribution. We then use this ball to run a naïve estimation procedure: clip the data to the ball $B_{1}$, then add noise to the empirical mean of the clipped data to obtain some initial estimate of the mean. Specifically, the noise will have magnitude proportional to $R_{1}/n$. Using this estimate, and knowledge of how we obtainesd it, we can draw a (hopefully significantly smaller) ball $B_{2}$ of radius $R_{2}$ that contains most of the mass and then repeat. After a few iterations, we will have some ball $B_{t}$ of radius $R_{t}$ that tightly contains most of the datapoints, and use this to make an accurate final private estimate of the mean with noise proportional to $R_{t}/n$. Our covariance estimation uses the same iterative approach, although the geometry is significantly more subtle.

2 Problem Formulation

An estimator $M(X)$ satisfies $\rho$-zCDP if for every pair of neighboring samples $X,X^{\prime}$ of size $n$, and every $\alpha\in(1,\infty)$, $D_{\alpha}(M(X)\|M(X^{\prime}))\leq\rho\alpha$, where $D_{\alpha}$ is the Rényi divergence of order $\alpha$.

This formulation sits in between general $(\varepsilon,\delta)$-differential privacy and the special case of $(\varepsilon,0)$-differential privacy, Formally, $(\varepsilon,0)$-DP $\Longrightarrow$ $\frac{1}{2}\varepsilon^{2}$-zCDP $\Longrightarrow$ $(\varepsilon\sqrt{2\log(1/\delta)}+\frac{1}{2}\varepsilon^{2},\delta)$-DP for every $\delta>0$ and better captures the privacy cost of private algorithms in high dimension [DSSU17].

We assume that our algorithms are given some a priori estimate of the mean in the form of a radius $R$ such that $\|\mu\|_{2}\leq R$ and some a priori estimate of the covariance in the form of $K$ such that $I\preceq\Sigma\preceq KI$, (equivalently all the singular values of $\Sigma$ lie between $1$ and $K$). We repeat that some a priori bound on $R,K$ is necessary for any algorithm that satisfies zCDP [BS16, KV18].

We measure the error in Mahalanobis distance $\|\cdot\|_{\Sigma}$, which reports how the error compares to the covariance of the distribution, and has the benefit of being invariant under affine transformations. Specifically,

For any distribution, the empirical mean $\hat{\mu}$ and empirical covariance $\hat{\Sigma}$ satisfy

and these estimators are minimax optimal. Our goal is to obtain estimators that have similar accuracy to the empirical mean and covariance. We note that the folklore naïve estimators (see e.g. [KV18, KLSU19]) for mean and covariance based on clipping the data to an appropriate ball and adding carefully calibrated noise to the empirical mean and covariance would guarantee

The main downside of the naïve estimators is their error increases rapidly with $R$ and $K$, and thus introduce large error unless the user has strong a priori knowledge of the mean and covariance. Requiring users to provide such a priori bounds is a major challenge in deployed systems for differentially private analysis (e.g. [GHK+16]). Our estimators have much better dependence on these parameters, both asymptotically and concretely.

For our theoretical analysis and most of our evaluation, we derive bounds on the error of our estimators assuming $P$ is specifically the Gaussian $N(\mu,\Sigma)$. Although, as we show in some of our experiments, our methods perform well even when we relax this assumption.

3 Related Work

The most relevant line of work is that initiated by Karwa and Vadhan [KV18], which studies private mean and variance estimation for Gaussian data, and focuses on important issues for practice like dealing with weak a priori bounds on the parameters. Later works studied the multivariate setting [KLSU19, CWZ19, KSSU19] and estimation under weaker moment assumptions [BS19, KSU20], though these investigations are primarily theoretical. Our algorithm for covariance estimation can be seen as a simpler and more practical variant of [KLSU19]. They provide an iterative procedure which, based on a privatized finds the subspaces of high and low variance they iteratively threshold eigenvalues to find directions of high and low variance, whereas we employ a softer method to avoid wasting information. One noteworthy work is [DFM+20], which provides practical private confidence intervals in the univariate setting. Instead, our investigation is focused on realizable algorithms for the multivariate setting. Several works consider private PCA or covariance estimation [DTTZ14, HP14, ADK+19], though, unlike our work, these methods assume strong a priori bounds on the covariance.

Other approaches for Gaussian estimation include [NRS07], which introduced the sample-and-aggregate paradigm, and [BKSW19] which employs a private hypothesis selection method. Zhang, Kamath, Kulkarni, and Wu privately estimate Markov Random Fields [ZKKW20], a generalization of product distributions over the hypercube. Dwork and Lei [DL09] introduced the propose-test-release framework for estimating robust statistics such as the median and interquartile range. For further coverage of private statistics, see [KU20].

Preliminaries

We begin by recalling the definition of differential privacy, and the variant of concentrated differential privacy that we use in this work.

A randomized algorithm $M:\mathcal{X}^{n}\rightarrow\mathcal{Y}$ satisfies $(\varepsilon,\delta)$-differential privacy ($(\varepsilon,\delta)$-DP) if for every pair of neighboring datasets $X,X^{\prime}\in\mathcal{X}^{n}$ (i.e., datasets that differ in exactly one entry),

When $\delta=0$, we say that $M$ satisfies $\varepsilon$-differential privacy or pure differential privacy.

A randomized algorithm $M:\mathcal{X}^{n}\rightarrow\mathcal{Y}$ satisfies $\rho$-zCDP if for every pair of neighboring datasets $X,X^{\prime}\in\mathcal{X}^{n}$,

where D_{\alpha}\mathopen{}\mathclose{{\left(M(X)||M(X^{\prime})}}\right) is the $\alpha$-Rényi divergence between $M(X)$ and $M(X^{\prime})$. Given two probability distributions $P,Q$ over $\Omega$, D_{\alpha}(P\|Q)=\frac{1}{\alpha-1}\log\mathopen{}\mathclose{{\left(\sum_{x}P(x)^{\alpha}Q(x)^{1-\alpha}}}\right).

Note that zCDP and DP are on different scales, but otherwise can be ordered from most-to-least restrictive. Specifically, $(\varepsilon,0)$-DP implies $\frac{\varepsilon^{2}}{2}$-zCDP, which implies roughly $(\varepsilon\sqrt{2\log(1/\delta)},\delta)$-DP for every $\delta>0$ [BS16].

Both these definitions are closed under post-processing and can be composed with graceful degradation of the privacy parameters.

If $M:\mathcal{X}^{n}\rightarrow\mathcal{Y}$ is $(\varepsilon,\delta)$-DP, and $P:\mathcal{Y}\rightarrow\mathcal{Z}$ is any randomized function, then the algorithm $P\circ M$ is $(\varepsilon,\delta)$-DP. Similarly if $M$ is $\rho$-zCDP then the algorithm $P\circ M$ is $\rho$-zCDP.

If $M$ is an adaptive composition of differentially private algorithms $M_{1},\dots,M_{T}$, then

if $M_{1},\dots,M_{T}$ are $(\varepsilon_{1},\delta_{1}),\dots,(\varepsilon_{T},\delta_{T})$-DP then $M$ is $(\sum_{t}\varepsilon_{t},\sum_{t}\delta_{t})$-DP, and

if $M_{1},\dots,M_{T}$ are $\rho_{1},\dots,\rho_{T}$-zCDP then $M$ is $(\sum_{t}\rho_{t})$-zCDP.

We can achieve differential privacy via noise addition proportional to sensitivity [DMNS06].

New Algorithms for Multivariate Gaussian Estimation

In this section, we present new algorithms for Gaussian parameter estimation. While these do not result in improved asymptotic sample complexity bounds, they will lead to algorithms which are much more practical in the multivariate setting. In particular, they will avoid the curse of dimensionality incurred by multivariate histograms, but also eschew many of the hyperparameters that arise in previous methods [KLSU19]. Note that our algorithms with $t=1$ precisely describe the naïve method that was informally outlined in Section 1.2. We describe the algorithms and sketch ideas behind the proofs, which appear in the appendix. Also in the appendix, we describe simpler univariate algorithms in the same style. Understanding these algorithms and proofs first might be helpful before approaching algorithms for the multivariate setting.

We first present our multivariate private mean estimation algorithm MVMRec (Algorithm 2). This is an iterative algorithm, which maintains a confidence ball that contains the true mean with high probability. For ease of presentation, we state the algorithm for a Gaussian with identity variance. However, by rescaling the data, the same argument works for an arbitrary known covariance $\Sigma$. In fact, the covariance doesn’t even need to be known exactly – we just need a proxy $\hat{\Sigma}$ such that $C_{1}I\preceq\hat{\Sigma}^{-1/2}\Sigma\hat{\Sigma}^{-1/2}\preceq C_{2}I$, where $0<C_{1}<C_{2}$ are absolute constants.

It remains to reason about MVM. We need to argue (a) privacy and (b) accuracy: given $B_{2}(c,r)\ni\mu$, it is likely to output $B_{2}(c^{\prime},r^{\prime})\ni\mu$; and c) progress: the radius $r^{\prime}$ output is much smaller than the $r$ input. The algorithm first chooses some $\gamma$ and clips the data to $B_{2}(c,r+\gamma)$. $\gamma$ is chosen based on Gaussian tail bounds such that if $\mu\in B_{2}(c,r)$, then none of the points will be affected by this operation. This bounds the sensitivity of the empirical mean, and applying the Gaussian mechanism guarantees privacy. While the noised mean will serve as a point estimate $c^{\prime}$, we actually have more: again using Gaussian tail bounds on the data combined with the added noise, we can define a radius $r^{\prime}$ such that $\mu\in B_{2}(c^{\prime},r^{\prime})$, establishing accuracy. Finally, a large enough $n$ will ensure $r^{\prime}<r/2$, establishing progress. Since each step reduces our radius $r_{i}$ by a constant factor, setting $t=O(\log R)$ will reduce the initial radius from $R$ to $O(1)$ as desired. Formalizing this gives the following theorem.

2 Multivariate Private Covariance Estimation

We describe our multivariate private covariance estimation algorithm MVCRec (Algorithm 4). The ideas are conceptually similar to mean estimation, but subtler due to the more complex geometry. We assume data is drawn from $N(0,\Sigma)$ where $I\preceq\Sigma\preceq KI$ for some known $K$. One can reduce to the zero-mean case by differencing pairs of samples. Further, if we know some PSD matrix $A$ such that $A\preceq\Sigma\preceq KA$ then we can rescale the data by $A^{-1/2}$.

To specify the algorithm we need a couple of tail bounds on the norm of points from a normal distribution, the spectral-error of the empirical covariance, and the spectrum of a certain random matrix. As these expressions are somewhat ugly, we define them outside of the pseudocode. In these expressions, we fix parameters $n,d,\beta_{s}$.

Similar to MVMRec, MVCRec repeatedly calls a private algorithm (MVC) that makes a constant-factor progress (with respect to some appropriate measure), and then runs the naïve algorithm (i.e., clip the data and noise the empirical covariance matrix). Rather than maintaining a ball containing the true mean, we maintain an ellipsoid (described via a PSD matrix) that upper bounds the true covariance in the Loewner order. For mathematical convenience, we work in a scaled version of the original space. That is, after each step, we rescale the problem so that this upper bound is the identity matrix, which simplifies reasoning about and describing the clipping procedure and noising mechanism. Progress holds with respect to the original problem in the unscaled domain: roughly speaking, either the upper bound on the variance in a direction decreases by a constant factor, or if this upper bound is already tight up to a constant factor, then the upper bound increases only slightly. As the upper and lower bounds on the variance in each direction are off by a factor of $K$, we show that $O(\log K)$ iterations suffice to get an upper bound which is at most a constant factor larger than the true covariance in each direction. At this point, we can apply the naïve clip-and-noise algorithm, which is accurate given enough samples.

The algorithm MVC is similar to MVM. We first clip the points at a distance based on Gaussian tail bounds with respect to the outer ellipsoid, which is unlikely to affect the dataset when it dominates the true covariance. After this operation, we can show that the sensitivity of an empirical covariance statistic is bounded using the following lemma. [KLSU19] proved a similar statement without an explicit constant, but the optimal constant is important in practice.

Applying the Gaussian mechanism (à la [DTTZ14]) in combination with this sensitivity bound, we again get a private point estimate $Z$ for the covariance, and can also derive a confidence ellipsoid upper bound. This time we require more sophisticated tools, including confidence intervals for the spectral norm of both a symmetric Gaussian matrix and the empirical covariance matrix of Gaussian data. Using a valid confidence ellipsoid ensures accuracy, and a sufficiently large $n$ again results in a constant factor squeezing of the ellipsoid, guaranteeing progress.

Putting together the analysis leads to the following theorem.

Conclusions

We provided the first effective and realizable algorithms for differentially private estimation of mean and covariance in the multivariate setting. We demonstrated that not only do these algorithms have strong theoretical guarantees (matching the state-of-the-art), but they are also accurate even at relatively low sample sizes and high dimensions. They significantly outperform all prior methods, possess few hyperparameters, and remain precise even when given minimal prior information about the data. In addition, we showed that our methods can be used for a private version of PCA, a task which is common in data science and exploratory data analysis. As we are seeing the rise of a number of new libraries for practical differentially private statistics and data analysis [The20, ABW20] we believe our results add an important tool to the toolkit for the multivariate setting.

Acknowledgments

GK thanks Aleksandar Nikolov for useful discussions about the proof of Lemma 3.2, and Argyris Mouzakis for pointing out a gap in a previous proof of Theorem 3.3.

References

Appendix A New Algorithms for Univariate Gaussian Parameter Estimation

In this section, we present our algorithms for estimating the mean and variance of a univariate Gaussian. We will write all our algorithms to give zCDP privacy guarantees, but the same approach can give pure DP algorithms in the univariate setting. One must simply switch Gaussian to Laplace noise, swap in the appropriate tail bounds for the confidence interval, and apply basic composition rather than zCDP composition.

We start with our univariate private mean estimation algorithm UVMRec (Algorithm 6). The guarantees are presented in Theorem A.1, note that the sample complexity is optimal in all parameters up to logarithmic factors [KV18]. While our algorithm and results are stated for a Gaussian with known variance, the same guarantees hold if the algorithm is only given the true variance up to a constant factor. Algorithm 6 is an iterative invocation of Algorithm 5, each step of which makes progress by shrinking our confidence interval for where the true mean lies. This is the simplest instantiation of our general algorithmic formula. Additionally, our proof of correctness is spelled out in full detail for this case – as the other proofs follow an almost identical structure, we only describe the differences.

We start by proving privacy. Observe that by application of the Gaussian mechanism (Lemma 2.6) in Line 4 of Algorithm 5, this algorithm is $\rho$-zCDP. Privacy of Algorithm 6 follows by composition of zCDP (Lemma 2.4).

We start by analyzing the $t-1$ iterations of the Line 3, each of which calls Algorithm 5. We prove two properties of this algorithm. Informally, it will always create a valid confidence interval, and the confidence interval shrinks by a constant factor. More formally:

A.2 Univariate Private Variance Estimation

Overall, the proof is very similar to that of Theorem A.1, so we only highlight the differences. First, we note that the proof of privacy is identical, via the Gaussian mechanism and composition of zCDP.

Appendix B Missing Proofs from Section 3

The proof is very similar to that of Theorem A.1, so we assume familiarity with that and only highlight the differences. First, we note that the proof of privacy is identical, via the Gaussian mechanism and composition of zCDP.

B.2 Proof of Lemma 3.2

Suppose the datasets differ in that one contains a point $X$ which is replaced by the point $Y$ in the other dataset.

The first inequality is since $\operatorname{Tr}(AB)\geq 0$, for any positive semi-definite matrices $A$ and $B$.

B.3 Proof of Theorem 3.3

Privacy again follows from the Gaussian mechanism and composition of zCDP. Note that this time, the sensitivity bound is not obvious – the analysis depends on Lemma 3.2.

To argue the utility guarantee of this procedure, we reason about the quantity $A_{i}$, which is the “scaling matrix” obtained at the end of the $i$th iteration. We rewrite $A_{i}$ as the product $M_{i}\times\dots\times M_{1}\times A_{0}=M_{i}A_{i-1}$, where $M_{i}$ is the matrix $U^{-1/2}$ obtained in Line 9 of the $i$th call to MVC. Specifically, we will argue that

for $i=0$ to $t-1$ with probability at least $1-\frac{\beta i}{2(t-1)}$. Note that, while the relationship $\Sigma\preceq(2-2^{-i})\Sigma+K2^{-i}I$ is trivial, the crucial aspects of this set of inequalities are that $A_{i}^{-1}A_{i}^{-1}$ is explicitly known after the $i$th call, serves as a valid upper bound for $\Sigma$, and is not too loose an upper bound on $\Sigma$.

We prove this by induction. Starting with $i=0$, we know that $A_{0}=\frac{1}{\sqrt{K}}I$ by definition, and thus $A_{0}^{-1}A_{0}^{-1}=KI$. By assumption in the theorem statement, we know that $\Sigma\preceq KI$. Furthermore, $KI\preceq\Sigma+KI$ trivially, and thus the base case holds with probability $1$.

Next, we take the inductive step, where we assume the statement holds for $A_{i-1}$, and prove it for $A_{i}$. By the inductive hypothesis, $A_{i-1}\Sigma A_{i-1}\preceq I$ with probability at least $1-\frac{\beta(i-1)}{2(t-1)}$, which we condition on. The analysis is similar to before: we bound the probability that a point gets adjusted in Line 3 using Fact C.2. We bound the spectral norm of the error due to sampling and the Gaussian noise matrix using Lemma C.3 and C.4, respectively. Combined, these give us that

with probability at least $1-\frac{\beta}{2(t-1)}$. Taking a union bound with the failure event from the induction hypothesis will give us the desired success probability of at least $1-\frac{\beta i}{2(t-1)}$, as the rest of the argument will be non-probabilistic in nature.

Focusing on the inequality $A_{i-1}\Sigma A_{i-1}\preceq U$, and multiplying both sides on the left and right by $A_{i-1}^{-1}$, we get

which is the first inequality we set out to prove.

It only remains to prove $A_{i}^{-1}A_{i}^{-1}\preceq(2-2^{-i})\Sigma+K2^{-i}I$. Given the expressions of $\eta,\nu$, and our choice of $n$, we can bound $\eta+\nu$ by $1/4$, and thus

We substitute this upper bound on $U$ into $A_{i}^{-1}A_{i}^{-1}$, giving

At this point, we apply the upper bound of the induction hypothesis:

Now, similar to before, we inspect the final call in Line 6. We have the following with probability at least $1-\beta/2$:

The first inequality is by assumption. The second and third are by (3) and the setting of $t=O(\log K)$. The final inequality uses the first inequality. Rearranging terms, we have that

With this in place, analysis follows similarly to Lemma 3.6 of [KLSU19]. Sketching the argument: our condition on $n$ implies that the Frobenius norm of both the empirical covariance and the noise added will be bounded by $\alpha$. Rescaling $Z$ by the scaling matrix $A_{t-1}$ gives the desired result.

Appendix C Concentration Inequalities and Tail Bounds

If $X\sim N(\mu,\sigma^{2})$, then $\Pr(|X-\mu|\geq\sigma\sqrt{2\log(2/\beta)})\leq\beta$.

If $X$ is a chi-squared random variable with $k$ degrees of freedom, then $\Pr(X-k\geq 2\sqrt{k\log(1/\beta)}+2\log(1/\beta))\leq\beta$ and $\Pr(k-X\geq 2\sqrt{k\log(1/\beta)})\leq\beta$. Thus, if $Y\sim N(0,I)$, then $\Pr(\|Y\|_{2}^{2}\geq d+2\sqrt{d\log(1/\beta)}+2\log(1/\beta))\leq\beta$.

We also need the following bound on the spectral error of an empirical covariance matrix.

Finally, we need a bounds on the spectral norm of a symmetric matrix with random Gaussian entries.

Let $Y$ be the $d\times d$ matrix where $Y_{ij}~N(0,\sigma^{2})$ for $i\leq j$, and $Y_{ij}=Y_{ji}$ for $i>j$. Then with probability at least $1-\beta$, we have the following bound:

First, we have the following bound on the expectation of the spectral norm:

This is from Theorem 1.1 of [BVH16], fixing the value of $\varepsilon=\frac{\log d}{d}$. The desired tail bound follows since the spectral norm is $2$-Lipschitz for this class of symmetric matrices, and by Gaussian concentration of Lipschitz functions (e.g., Proposition 5.34 of [Ver12]). ∎