Differentially Private Bayesian Inference for Exponential Families

The study of private inference has been sparked by growing concern regarding the analysis of data when it stems from sensitive sources. We present the first method for private Bayesian inference in exponential families that properly accounts for noise introduced by the privacy mechanism. It is efficient because it works only with sufficient statistics and not individual data. Unlike other methods, it gives properly calibrated posterior beliefs in the non-asymptotic data regime.


Introduction
Differential privacy is the dominant standard for privacy [1]. A randomized algorithm that satisfies differential privacy offers protection to individuals by guaranteeing that its output is insensitive to changes caused by the data of any single individual entering or leaving the data set. An algorithm can be made differentially private by applying one of several general-purpose mechanisms to randomize the computation in an appropriate way, for example, by adding noise calibrated to the sensitivity of the quantity being computed, where sensitivity captures how much the quantity depends on any individual's data [1]. Due to the obvious importance of protecting individual privacy while drawing population level inferences from data, differentially private algorithms have been developed for a broad range of machine learning tasks [2][3][4][5][6][7][8][9].
There is a growing interest in private methods for Bayesian inference [10][11][12][13][14]. In Bayesian inference, a modeler selects a prior distribution p(θ) over some parameter, observes data x that depends probabilistically on θ through a model p(x | θ), and then reasons about θ through the posterior distribution p(θ | x), which quantifies updated beliefs and uncertainty about θ after observing x. Bayesian inference is a core machine learning task and there is an obvious need to be able to conduct it in a way that protects privacy when x is sensitive. Additionally, recent work has identified surprising connections between sampling from posterior distributions and differential privacy-for example, a single perfect sample from p(θ | x) satisfies differential privacy for some setting of the privacy parameter [10][11][12][13].
An "obvious" way to conduct private Bayesian inference is to privatize the computation of the posterior, that is, to design a differentially private algorithm A that outputs y = A(x) with the goal that y ≈ p(θ | x) is a privatized representation of the posterior. However, using y directly as "the 32nd Conference on Neural Information Processing Systems (NIPS 2018), Montréal, Canada. posterior" will not correctly quantify beliefs, because the Bayesian modeler never observes x, she observes y; her posterior beliefs are now quantified by p(θ | y). This paper will take a different approach to private Bayesian inference by designing a pair of algorithms: The release mechanism A computes a private statistic y = A(x) of the input data; the inference algorithm P computes p(θ | y). These algorithms should satisfy the following criteria: • Privacy. The release mechanism A is differentially private. By the post-processing property of differential privacy [15], all further computations are also private. • Calibration. The inference algorithm P can efficiently compute or approximate the correct posterior, p(θ | y) (see Section 4 for our process to measure calibration). • Utility. Informally, the statistic y should capture "as much information as possible" about x so that p(θ | y) is "close" to p(θ | x).
Importantly, the release mechanism A is public, so the distribution p(y | x) is known. Williams and McSherry first suggested conducting inference on the output of a differentially private algorithm and showed how to do this for the factored exponential mechanism [16]; see also [17][18][19][20].
Our work focuses specifically on Bayesian inference when the private data X = x 1:n is an iid sample of (publicly known) size n from an exponential family model p(x i | θ). Exponential families include many of the most familiar parametric probability models. We will adopt a straightforward release mechanism where the Laplace mechanism [1] is used to release noisy sufficient statistics y [12,19], which are a finite-dimensional quantity that capture all the information about θ [21].
The technical challenge is then to develop an efficient general-purpose inference algorithm P. One challenge is computational efficiency. The exact posterior p(θ | y) ∝ p(θ)p(x 1:n | θ)p(y|x 1:n )dx 1:n integrates over all possible data sets [16], which is intractable to do directly for large n. We integrate instead over the sufficient statistics s, which have fixed dimension and completely characterize the posterior; furthermore, since they are a sum over individuals, p(s | θ) is asymptotically normal. We develop an efficient Gibbs sampler that uses a normal approximation for s together with variable augmentation to model the Laplace noise in a way that yields simple updates [22].
A second challenge is that the sufficient statistics may be unbounded, which makes their release incompatible with the Laplace mechanism. We address this by imposing truncation bounds and only computing statistics from data that fall within the bounds. We show how to use automatic differentiation and a "random sum" central limit theorem to compute the parameters of the normal approximation p(s | θ) for a truncated exponential family when the number of individuals that fall within the truncation bounds is unknown.
Our overall contribution is the pairing of an existing simple release mechanism A with a novel, efficient, and general-purpose Gibbs sampler P that meets the criteria outlined above for private Bayesian inference in any univariate exponential family or multivariate exponential family with bounded sufficient statistics. 1 We show empirically that when compared with competing methods, ours is the only one that provides properly calibrated beliefs about θ in the non-asymptotic regime, and that it provides good utility compared with other private Bayesian inference approaches.

Differential Privacy
Differential privacy requires that an individual's data has a limited effect on the algorithm's behavior. In our setting, a data set X = x 1:n := (x 1 , . . . , x n ) consists of records from n individuals, where x i ∈ R d is the data of the ith individual. We will assume n is known. Differential privacy reasons about the hypothesis that one individual chooses not to remove their data from the data set, and their record is replaced by another one. 2 Let nbrs(X) denote the set of data sets that differ from X by exactly one record-i.e., if X ∈ nbrs(X), then X = (x 1:i , x i , x i+1:n ) for some i. 1 There are remaining technical challenges for multivariate models with unbounded sufficient statistics that we leave for future work. 2 This variant assumes n remains fixed, which is sometimes called bounded differential privacy [23].
Definition 1 (Differential Privacy; Dwork et al. [1]). A randomized algorithm A satisfies -differential privacy if for any input X, any X ∈ nbrs(X) and any subset of outputs O ⊆ Range(A), We achieve differential privacy by injecting noise into statistics that are computed on the data. Let f be any function that maps datasets to R d . The amount of noise depends on the sensitivity of f .
We drop the subscript f when it is clear from context. Our approach achieves differential privacy through the application of the Laplace mechanism. Definition 3 (Laplace Mechanism; Dwork et al. [1]). Given a function f that maps data sets to R m , the Laplace mechanism outputs the random variable L(X) ∼ Lap f (X), ∆ f / from the Laplace distribution, which has density Lap(z; u, b) = (2b) −m exp (− z − u 1 /b). This corresponds to adding zero-mean independent noise u i ∼ Lap(0, ∆ f / ) to each component of f (X).
A final important property of differential privacy is post-processing [15]; if an algorithm A isdifferentially private, then any algorithm that takes as input only the output of A, and does not use the original data set X, is also -differentially private.

Private Bayesian Inference in Exponential Families
We consider the canonical setting of Bayesian inference in an exponential family. The modeler posits a prior distribution p(θ), assumes the data x 1:n is an iid sample from an exponential family model p(x | θ), and wishes to compute the posterior p(θ | x 1:n ). An exponential family in natural parameterization has density where η are the natural parameters, t(x) is the sufficient statistic, A(η) = h(x) exp η T t(x) dx is the log-partition function, and h(x) is the base measure. The density of the full data is where h(x 1:n ) = n i=1 h(x i ) and t(x 1:n ) = n i=1 t(x i ). Notice that once normalizing constants are dropped, this density is dependent on the data only directly through the sufficient statistics, s = t(x 1:n ).
We will write exponential families more generally as p(x | θ) to indicate the case when the natural parameters η = η(θ) depend on a different parameter vector θ.
Every exponential family distribution has a conjugate prior distribution p(θ; λ) [24] with hyperparameters λ. A conjugate prior has the property that, if it is used as the prior, then the posterior belongs to the same family, i.e., p(θ | x 1:n ; λ) = p(θ; λ ) for some λ that depends only on λ, n, and the sufficient statistics s. We write this function as λ = Conjugate-Update(λ, s, n); our methods are not tied to the specific choice of conjugate prior, only that the posterior parameters can be calculated in this form. See supplementary material for a general form of Conjugate-Update.

Release Algorithm: Noisy Sufficient Statistics
If privacy were not a concern, the Bayesian modeler would simply compute the sufficient statistics s = t(x 1:n ) and use them to update the posterior beliefs. However, to maintain privacy, the modeler must access the sensitive data only through a randomized release mechanism A. As a result, in order to obtain proper posterior beliefs the modeler must account for the randomization of the release mechanism by performing inference.
We take the simple approach of releasing noisy sufficient statistics via the Laplace mechanism, as in [12,13,19]. Sufficient statistics are a natural quantity to release. They are an "information bottleneck"-a finite-dimensional quantity that captures all the relevant information about θ. The released value is y = A( When t(·) is unbounded this quantity becomes infinite; we will modify the release mechanism so the sensitivity is finite (Sec. 3.3).

Basic Inference Approach: Bounded Sufficient Statistics
The goal of the inference algorithm P is to compute p(θ | y). We first develop the basic approach for the simpler case when t(x) is bounded, and then extend both A and P to handle the unbounded case. The full joint distribution of the probability model can be expressed as: where p(θ) = p(θ; λ) is a conjugate prior and the goal is to compute a representation of p(θ | y) ∝ s p(θ, s, y)ds by integrating over the sufficient statistics. We will develop a Gibbs sampler to sample from this distribution. There are two main challenges. First, the distribution p(s | θ) is obtained by marginalizing over the data sample x 1:n , and is usually not known in closed form. We will address this with an asymptotically correct normal approximation. Second, when resampling s within the Gibbs algorithm, we require the full conditional distribution of s given the other variables, which is proportional to p(s|θ)p(y | s). Care must be taken to make it easy to sample from this conditional distribution. We address this via variable augmentation. We discuss our approach to both challenges in detail below.
Normal approximation of p(s | θ). The exact form of the sufficient statistic distribution p(s | θ) is obtained by marginalizing over the data: In general, the exact form of this distribution is not available. In some cases, it is-for example if x ∼ Bernoulli(θ) then s ∼ Binomial(n, θ)-but even then it may not lead to a tractable full conditional for s.
Properties of exponential families pave the way toward a general approach that always leads to a tractable full conditional. By the central limit theorem (CLT), because s = i t(x i ) is a sum of iid random variables, it is asymptotically normal. It can be approximated as are the mean and variance of the sufficient statistic of a single individual. This approximation is asymptotically correct [25]. The quantities µ and Σ can be computed using well-known properties of exponential families [25]: where η = η(θ) is the natural parameter.
Note that we will not use this approximation for Gibbs updates of θ. Instead, we will compute the conditional p(θ | s) using standard conjugacy formulas. In this sense, we maintain two views of the joint distribution p(θ, s)-when updating θ, it is the standard exponential family model, which leads to conjugate updates; when updating s, it is approximated as p(θ)N (s; nµ, sΣ), which will lead to simple updates when combined with a variable augmentation technique.
Variable augmentation for p(y | s). We seek a tractable form for the full conditional of s under the normal approximation, which is the product of a normal density and a Laplace density: A similar situation arises in the Bayesian Lasso [22], and we will employ the same variable augmentation trick. A Laplace random variable z ∼ Lap(u, b) can be written as a scale mixture of normals by introducing a latent variable σ 2 ∼ Exp(1/(2b 2 )), i.e., the distribution with density 1/(2b 2 ) exp −σ 2 /(2b 2 ) and letting z ∼ N (u, σ 2 ). We apply this separately to each dimension of the vector y so that: s ∼ NormProduct nµ, nΣ, y, diag(σ 2 ) 6: The Gibbs Sampler. After the normal approximation and variable augmentation, the generative process is as shown to the right. The final Gibbs sampling algorithm is shown in Algorithm 1. Note that the update for θ is based on conjugacy in the exact distribution p(θ, s), while the update for s uses the density of the generative process to the right, so that p(s | θ, σ 2 , y) ∝ p(s | θ) p(y | σ 2 , s), which is a product of two normal densities where µ s and Σ s are are defined in Algorithm 1 [26]. The update for σ 2 follows Park and Casella [22]; the inverse Gaussian density is InverseGaussian( . Full derivations are given in the supplement.

Unbounded Sufficient Statistics and Truncated Exponential Families
The Laplace mechanism does not apply when the sufficient statistics are unbounded, because ∆ s = max x,y t(x) − t(y) 1 = ∞. Thus, we need a new release mechanism A and inference algorithm P.
We present a solution for the case when x is univariate. All elements of the solution can generalize to higher dimensions, except that one step will have running time that is exponential in d; we leave improvement of this to future work and focus on the simpler univariate case.
Release mechanism. Our solution is to truncate the support of the (now univariate) p(x | θ) to x ∈ [a, b], where a and b are finite bounds provided by the modeler. If the modeler cannot select bounds a priori, they may be selected privately as a preliminary step using a variant of the exponential mechanism (see PrivateQuantile in Smith [27]). 3 Then, given truncation bounds, the data owner redacts individuals where x i / ∈ [a, b] and reports the truncated sufficient statisticŝ is the indicator function of the set S. The sensitivity ofŝ is now ∆ŝ = max x,y∈R t (x) −t(y) 1 wheret(x) = 1 [a,b] (x) t(x). An easy upper bound for this quantity (see supplement) is: where t j (x) is the jth component of the sufficient statistics. The bounds [a, b] will be selected so this quantity is bounded. The released value is y ∼ Lap(ŝ, ∆ŝ/ ).
Inference: Truncated Exponential Family. Several new challenges arise for inference. The quantityŝ is no longer a sufficient statistic for the model p(x | θ), and we will need new insights to understand p(ŝ | θ) and p(θ |ŝ). Sinceŝ is a sum over individuals where x i ∈ [a, b], it will be useful to examine the probability of the event x ∈ [a, b] as well as the conditional distribution of x given this event. To facilitate a general development, assume a generic truncation interval [v, w], not necessarily equal to [a, b]. Let F (x; θ) = x −∞ p(x | θ)dx be the CDF of the original (univariate) exponential family model. It is clear that Pr(x ∈ [v, w]) = F (w; θ) − F (v; θ). The conditional distribution of x given x ∈ [v, w] is a truncated exponential family, which, in its natural parameterization is: Note that this is still an exponential family model (with a modified base measure), and all of the standard results apply, such as the existence of a conjugate prior and the formulas in Eq. (1) for the mean and variance of t(x) under the truncated distribution.
Random sum CLT for p(ŝ | θ). We would like to again apply an asympotic normal approximation forŝ, but we do not know how many individuals fall within the truncation bounds. The "random sum CLT" of Robbins [28] applies to the setting where the number of terms in the sum is itself a random variable. The sum can be rewritten asŝ = Proof. Each term of the sum has meanμ and varianceΣ, and the number of terms is N ∼ Binom(n, q). The result follows from Robbins [28].
Computingμ andΣ by automatic differentiation (autodiff). To use the normal approximation we need to computeμ andΣ. Lemma 1. Let p(x | θ) be a univariate exponential family model and letp(x | θ) be the corresponding exponential family model truncated to generic interval [v, w]. Then Proof. It is straightforward to derive from Eq. (2) thatÂ(η) = A(η) + log F (w; η) − F (v; η) . The result follows from applying Eq. (1) to this expression forÂ(η). See the supplement for derivation of A(η) and proof of this lemma.
We will use Equations (3) and (4)  When x is multivariate, analogous expressions can be derived forμ andΣ. The adjustment factors will include multivariate CDFs, with a number of terms that grow exponentially in d. This is currently the main limitation in applying our methods to multivariate models with unbounded sufficient statistics.
Conjugate updates for p(θ |ŝ). The final issue is the distribution p(θ |ŝ), which is no longer characterized by conjugacy becauseŝ are not the full sufficient statistics. We again turn to variable Conditioned on all other variables, each component is multivariate normal, so the sum s is also multivariate normal. We can therefore sample s and then sample from p(θ | s) using conjugacy. We will also need to drawŝ c separately to be used to update σ 2 .
The Gibbs Sampler. The (approximate) generative process in the unbounded case is: The Gibbs sampler to sample from this distribution is given in Algorithm 2. Note that in Line 8 we employ rejection sampling in which sufficient statistics are sampled until the values drawn are valid for the given data model, e.g., s must be positive for the binomial distribution. The RS-CLT algorithm to compute parameters of the random sum CLT is shown in Algorithm 3.

Experiments
We design experiments to measure the calibration and utility of our method for posterior inference.
We conduct experiments for the binomial model with beta prior, the multinomial model with Dirichlet prior, and the exponential model with gamma prior. The last model is unbounded and requires truncation; we set the bounds to keep the middle 95% of individuals, which is reasonable to assume known a priori for some cases, such as modeling human height.
Methods. We run our Gibbs sampler for 5000 iterations after 2000 burnin iterations (see supplementary material for convergence results), which we compare to two baselines. The first method uses the same release mechanism as our Gibbs sampler and performs conjugate updates using the noisy sufficient statistics [12,13]. This method converges to the true posterior as n → ∞ because the Laplace noise will eventually become negligible compared to sampling variability [12]. However, the noise is not negligible for moderate n; we refer to this method as "naive". For truncated models we allow the naive method to "cheat" by accessing the noisy untruncated sufficient statistics s. Thus the method is not private, and receives strictly more information than our Gibbs sampler, but with the same magnitude noise. This allows us to demonstrate miscalibration without highly technical modifications to the baseline method to be able to deal with truncated sufficient statistics.
The second baseline is a version of the one-posterior sampling (OPS) mechanism [11][12][13], which employs the exponential mechanism [29] to release samples from a privatized posterior. We release 100 samples using the method of [12], each with ops = /100, such that the entire algorithm achieves -differential privacy. Private MCMC sampling [11] is a more sophisticated method to release multiple samples from a privatized posterior and could potentially make better use of the privacy budget; however, private MCMC will also necessarily be miscalibrated, and only achieves the weaker privacy guarantee of ( , δ)-differential privacy for δ > 0, so would not be direct comparable to our method. OPS serves as a suitable baseline that achieves -differential privacy. We include OPS only for experiments on the binomial model, for which it requires the support of θ to be truncated to [a 0 , 1 − a 0 ] where a 0 > 0. We set a 0 = 0.1.
We also include a non-private posterior for comparison, which performs conjugate updates using the non-noisy sufficient statistics.
Evaluation. We evaluate both the calibration and utility of the posterior. For calibration we adapt a method of Cook et al. [30]: the idea is to draw iid samples (θ i , x i ) from the joint model p(θ)p(x | θ), and conduct posterior inference in each trial. Let F i (θ) be the CDF of the true posterior p(θ | x i ) in trial i. Then we know that In other words, the actual parameter θ i is equally likely to land at any quantile of the posterior. To test the posterior inference procedure, we instead compute U i as the quantile at which θ i lands within a set of samples from the approximate posterior. After M trials of the whole procedure we test for uniformity of U 1:M using the Kolmogorov-Smirnov goodnessof-fit test [31], which measures the maximum distance between the empirical CDF of U 1:M and the uniform CDF; lower values are better and zero corresponds to perfect uniformity. We also visualize the empirical CDFs to assess calibration qualitatively.
Higher utility of a private posterior is indicated by closeness to the non-private posterior, which we measure with maximum mean discrepancy (MMD), a kernel-based statistical test to determine if two sets of samples are drawn from different distributions [32]. Given m i.i.d. samples (p, q) ∼ P × Q, an unbiased estimate of the MMD is where k is a continuous kernel function; we use a standard normal kernel. The higher the value the more likely the two samples are drawn from different distributions.
Results. Figure 1a shows the results for three models and varying n and . Our method (Gibbs) achieves the same calibration level as non-private posterior inference for all settings. The naive method ignores noise and is too confident about parameter values implied by treating the noisy sufficient statistics as true ones; it is only well-calibrated with increasing n and when noise becomes negligible relative to population size. OPS is not calibrated because it samples from an over-dispersed version of p(θ | x). Figure 1b shows the empirical CDF plots for n = 1000 and = 0.01. Our method and the non-private method are both perfectly calibrated. The naive method's over-confidence in the wrong sufficient statistics causes its posterior to usually be too tight at the wrong value; thus the true parameter always lies in a tail of the approximate posterior, so too much mass is placed near 0 and 1. OPS shows the opposite behavior: its posterior is always too diffuse, so the true parameter lies close to the middle. For multinomial we show measures only for the parameter of the first category, but results hold for all categories. Figure 1c shows the MMD test statistic between each method and the non-private posterior, used as a measure of utility. Our method consistently achieves utility at least as good as the naive method for binomial and multinomial models. We omit OPS, which is never calibrated. For the exponential model (not shown) we did not obtain conclusive utility comparisons due to the lack of a naive baseline that properly handles truncation; the "cheating" naive method from our calibration experiments sometimes attains higher utility than our method, and sometimes lower, but this comparison is not meaningful because it receives strictly more information.

Conclusion
We presented a Gibbs sampling approach for private posterior inference in exponential family models. Rather than trying to approximate the posterior of p(θ | x 1:n ), we divide our procedure into a private release mechanism y = A(x 1:n ) and an inference algorithm P that computes p(θ | y). The release mechanism is designed to facilitate inference. We develop a general-purpose Gibbs sampler that applies to any exponential family model that has bounded sufficient statistics; a truncated version applies to univariate models with unbounded sufficient statistics. The Gibbs sampler uses general properties of exponential families to approximate the distribution of the sufficient statistics, and therefore avoids the need to reason about individuals. Promising lines of future work are to develop efficient methods for multivariate exponential families with unbounded sufficient statistics, and to develop methods for conditional models based on exponential families, such as generalized linear models.

A Properties of Exponential Families
A.1 Form of Conjugate-Update(λ, x 1:n ) Following Diaconis and Ylvisaker [24], the prior is where the parameters are λ = [λ 1 , λ 2 ] and sufficient statistics are [η, −A(η)] The posterior after observing x 1:n is Proof: The proof for Varp[t(x)] is similar.

B Derivation of σ 2 Gibbs update
We fully derive the Gibbs update for the noise variance σ 2 of the augmented model as stated in Park and Casella [22]. We represent the Laplace distribution with scale b = ∆ s / as a scale mixture of normals, i.e. a zero-mean normal with an exponential prior on the variance:

C Sensitivity of Sufficient Statistics in Truncated Model
Recall thatt(x) = 1 [v,w]  |t j (x)|, max x,y∈ [v,w] |t j (x) − t j (y)| D Proof of uniformity of CDF transform used by Cook et al. [30] Claim: Let X be a random variable with CDF F . The random variable U = F (X) is uniformly distributed.

Proof:
Pr(U ≤ u) = Pr(F (X) ≤ u) = Pr F −1 (F (X)) ≤ F −1 (u) = Pr(X ≤ F −1 (u)) = F F −1 (u) = u E Convergence of Gibbs Sampler Figure 2 shows the progress of sampled model parameters over the course of 500 iterations for both binomial and exponential models. For both models the samples quickly converge to the vicinity of the true parameter.