# Differentially Private Release of High-Dimensional Datasets using the Gaussian Copula

←

**Page content transcription**

If your browser does not render page correctly, please read the page content below

Differentially Private Release of High-Dimensional Datasets using the Gaussian Copula Hassan Jameel Asghar1,2 , Ming Ding1 , Thierry Rakotoarivelo1 , Sirine Mrabet1 , Mohamed Ali Kaafar1,2 1 Data61, CSIRO, Sydney, Australia arXiv:1902.01499v1 [cs.CR] 4 Feb 2019 {hassan.asghar, ming.ding, thierry.rakotoarivelo, sirine.mrabet, dali.kaafar}@data61.csiro.au 2 Macquarie University, Sydney, Australia {hassan.asghar, dali.kaafar}@mq.edu.au February 6, 2019 Abstract We propose a generic mechanism to efficiently release differentially private synthetic versions of high- dimensional datasets with high utility. The core technique in our mechanism is the use of copulas, which are functions representing dependencies among random variables with a multivariate distribution. Specif- ically, we use the Gaussian copula to define dependencies of attributes in the input dataset, whose rows are modelled as samples from an unknown multivariate distribution, and then sample synthetic records through this copula. Despite the inherently numerical nature of Gaussian correlations we construct a method that is applicable to both numerical and categorical attributes alike. Our mechanism is efficient in that it only takes time proportional to the square of the number of attributes in the dataset. We propose a differentially private way of constructing the Gaussian copula without compromising computa- tional efficiency. Through experiments on three real-world datasets, we show that we can obtain highly accurate answers to the set of all one-way marginal, and two-and three-way positive conjunction queries, with 99% of the query answers having absolute (fractional) error rates between 0.01 to 3%. Furthermore, for a majority of two-way and three-way queries, we outperform independent noise addition through the well-known Laplace mechanism. In terms of computational time we demonstrate that our mechanism can output synthetic datasets in around 6 minutes 47 seconds on average with an input dataset of about 200 binary attributes and more than 32,000 rows, and about 2 hours 30 mins to execute a much larger dataset of about 700 binary attributes and more than 5 million rows. To further demonstrate scalability, we ran the mechanism on larger (artificial) datasets with 1,000 and 2,000 binary attributes (and 5 mil- lion rows) obtaining synthetic outputs in approximately 6 and 19 hours, respectively. These are highly feasible times for synthetic datasets, which are one-off releases. Keywords. differential privacy, synthetic data, high dimensional, copula 1 Introduction There is an ever increasing demand to release and share datasets owing to its potential benefits over controlled access. For instance, once data is released, data custodian(s) need not worry about access controls and continual support. From a usability perspective, data release is more convenient for users (expert analysts and novices alike) as compared to access through a restricted interface. Despite its appeal, sharing datasets, especially when they contain sensitive information about individuals, has privacy implications which have been well documented. Current practice, therefore, suggests privacy-preserving release of datasets. Ad hoc techniques such as de-identification, which mainly rely on properties of datasets and assumptions on what background information is available, have failed to guarantee privacy Narayanan and Shmatikov [2008], Ohm [2009]. Part of the reason for the failure is the lack of a robust definition of privacy underpinning these techniques. 1

This gave rise to the definition of differential privacy Dwork et al. [2006]. Differential privacy ties the privacy property to the process or algorithm (instead of the dataset) and, informally, requires that any output of the algorithm be almost equally likely even if any individual’s data is added or removed from the dataset. A series of algorithms have since been proposed to release differentially private datasets, often termed as synthetic datasets. Simultaneously, there are results indicating that producing synthetic datasets which accurately answer a large number of queries is computationally infeasible, i.e., taking time exponential in the number of attributes (dimension of the dataset) Ullman and Vadhan [2011]. This is a serious roadblock as real-world datasets are often high-dimensional. However, infeasibility results in Ullman and Vadhan [2011], Ullman [2013] are generic, targeting provable utility for any input data distribution. It may well be the case that a large number of real-world datasets follow constrained distributions which would make it computationally easier to output differentially private synthetic datasets that can accurately answer a larger number of queries. Several recent works indicate the plausibility of this approach claiming good utility in practice Zhang et al. [2014], Li et al. [2014]. Ours is a continuation of this line of work. In this paper, we present a generic mechanism to efficiently generate differentially private synthetic versions of high dimensional datasets with good utility. By efficiency, we mean that our mechanism can output a synthetic dataset in time O(m2 n), where m is the total number of attributes in the dataset and n the total number of rows. Recall that the impossibility results Ullman and Vadhan [2011] suggest that algorithms for accurately answering a large number of queries are (roughly) expected to run in time poly(2m , n). Thus, our method is scalable for high dimensional datasets, i.e., having a large m. In terms of utility, our generated synthetic dataset is designed to give well approximated answers to all one-way marginal and two-way positive conjunction queries. One-way marginal queries return the number of rows in the dataset exhibiting a given value x or its negation x (not x) under any attribute X. Similarly, two-way positive conjunction queries return the number of rows that exhibit any pair of values (x, y) under an attribute pair (X, Y ). This forms a subset of all two-way margins; the full set also includes negations of values, e.g., rows that satisfy (x, y). While this may seem like a small subset of queries, there are results showing that even algorithms that generate synthetic datasets to (accurately) answer all two-way marginal queries are expected to be computationally inefficient Ullman and Vadhan [2011], Ullman [2013]. Furthermore, we show that our mechanism provides good utility for other queries as well by evaluating the answers to 3-way positive conjunction queries on the synthetic output. The key to our method is the use of copulas Nelsen [2006] to generate synthetic datasets. Informally, a copula is a function that maps the marginal distributions to the joint distribution of a multivariate distribu- tion, thus defining dependence among random variables. Modelling the rows of the input dataset as samples of the (unknown) multivariate distribution, we can use copulas to define dependence among attributes (modelled as univariate random variables), and finally use it to sample rows from the target distribution and generate synthetic datasets. Specifically, we use the Gaussian copula [Nelsen, 2006, p. 23], which defines the dependence between one-way margins through a covariance matrix. The underlying assumption is that the relationship between different attributes of the input dataset is completely characterised by pairwise covariances. While this may not be true in practice, it still preserves the correlation of highly correlated attributes in the synthetic output. Our main reason for using the Gaussian copula is its efficiency, as its run-time is proportional to square of the number of attributes. We remark that our use of the Gaussian copula does not require the data attributes to follow a Gaussian distribution. Only their dependencies are assumed to be captured by the Gaussian copula. Importantly, as claimed, our method is generic. This is important since an input dataset is expected to be a mixture of numerical (ordinal) and categorical (non-ordinal) attributes meaning that we cannot use a unified measure of correlation between attributes. A common technique to work around this is to create an artificial order on the categorical attributes Iyengar [2002]. However, as we show later, the resulting correlations are inherently artificial and exhibit drastically different results if a new arbitrary order is induced. Our approach is to convert the dataset into its binary format (see Section 2.2) and then use a single correlation measure, Pearson’s product-moment correlation, for all binary attributes. This eliminates the need for creating an artificial order on the categorical attributes. Our method is thus more generic than a similar method proposed in Li et al. [2014], which only handles attributes with large (discrete) domains and induces artificial order on categorical attributes.1 1 See Section 5 for a further discussion on the differences between the two works. 2

We experimentally evaluate our method on three real-world datasets: the publicly available Adult dataset Lichman [2013] containing US census information, a subset of the social security payments dataset provided to us by the Department of Social Services (DSS), a department of the Government of Australia,2 , and a hospital ratings dataset extracted from a national patient survey in the United States, which we call the Hospital dataset. The Adult dataset consists of 14 attributes (194 in binary) and more than 32,000 rows, the DSS dataset has 27 attributes (674 in binary) and more than 5,000,000 rows, and the Hospital dataset contains 9 attributes (1,201 in binary) and 10,000 rows. Generation of synthetic datasets took around 6 minutes 47 seconds (on average) for the Adult dataset, around two and a half hours on average for the DSS dataset, and 46 minutes on average for the Hospital dataset. To further check the scalability of our method, we ran it on two artificial datasets each with more than 5 million rows and 1,000 and 2,000 binary attributes, resulting in run-times of approximately 6 and 19 hours, respectively. Since synthetic data release is a one-off endeavour, these times are highly feasible. In terms of utility, we show that 99% of the queries have an absolute error of less than 500 on Adult, less than 1,000 on DSS and around 300 on the Hospital dataset, where absolute error is defined as the absolute difference between true and differentially private answers computed from the synthetic dataset. In terms of two-way positive conjunction queries, we again see that 99% of the queries have an absolute error of around 500 for Adult, less than 500 for DSS and only around 20 for the Hospital dataset. Furthermore, for most of the two-way queries we considerably outperform the Laplace mechanism Dwork et al. [2006] of adding independent Laplace noise to each of the queries. Note that a further advantage of our method is that unlike the Laplace mechanism, we generate a synthetic dataset. We further expand our utility analysis to include three-way queries and show that our method again calibrates noise considerably better than the aforementioned Laplace mechanism with 99% of the queries having absolute error of around 400 for both Adult and DSS, and only around 50 for the Hospital dataset. Our utility analysis is thorough; we factor in the possibility that real-world datasets may have significantly high number of uncorrelated attributes which implies that a synthetic data generation algorithm might produce accurate answers by chance (for the case of two-way or higher order marginals). Thus we separate results for highly correlated and uncorrelated attributes to show the accuracy of our mechanism. Perhaps one drawback of our work is the lack of a theoretical accuracy bound; but this is in line with many recent works Zhang et al. [2014], Cormode et al. [2012], Li et al. [2014] which, like us, promise utility in practice. The rest of the paper is organized as follows. We give a brief background on differential privacy and copulas in Section 2. We describe our mechanism in Section 3. Section 4 contains our experimental utility and performance analysis. We present related work in Section 5, and conclude in Section 6. 2 Background concepts 2.1 Notations We denote the original dataset by D, which is modelled as a multiset of n rows from the domain X = A1 × · · · × Am , where each Aj represents an attribute, for a total of m attributes. We assume the number of rows n is publicly known. We denote the set of all n-row databases as X n . Thus, D ∈ X n . The ith row of D (i) (i) (i) (i) is denoted as (X1 , X2 , . . . , Xm ), where Xj ∈ Aj . In the sequel, where discussing a generic row, we will drop the superscript for brevity. The notation R represents the real number line [−∞, ∞], and I denotes the interval of real numbers [0, 1]. The indicator function I{P } evaluates to 1 if the predicate P is true, and 0 otherwise. 2.2 Dummy Coding A key feature of our method is to convert the original dataset into its binary format through dummy coding, where each attribute value of Aj ∈ D is assigned a new binary variable. For instance, consider the “country” attribute having attribute values USA, France and Australia. Converted to binary we will have three columns (binary attributes) one for each of the three values. For any row, a value of 1 corresponding to any of these binary columns indicates that the individual is from that particular country. Thus, we have a 2 https://www.dss.gov.au/. 3

Pm total of d = i=1 |Ai | binary attributes in the binary version DB of the dataset D. In case of continuous valued attributes, this is done by first binning the values into discrete bins. Thus for a continuous valued (i) (i) (i) attribute A, |A| indicates the number of bins. The ith row in DB is denoted as (X1 , X2 , . . . , Xd ). Note 3 that while the mapping from D to DB is unique, the converse is not always true. Note further that this way of representing datasets is also known as histogram representation [Dwork and Roth, 2014, §2.3]. 2.3 Overview of Differential Privacy Two databases D, D0 ∈ X n are called neighboring databases if they differ in only one row. Definition 1 (Differential privacy Dwork et al. [2006], Dwork and Roth [2014]). A randomized algorithm (mechanism) M : X n → R is (, δ)-differentially private if for every S ⊆ R, and for all neighbouring databases D, D0 ∈ X n , the following holds P(M(D) ∈ S) ≤ e P(M(D0 ) ∈ S) + δ. If δ = 0, then M is -differentially private. The mechanism M might also take some auxiliary inputs such as a query or a set of queries (to be defined shortly). The parameter δ is required to be a negligible function of n [Dwork and Roth, 2014, §2.3, p. 18],[Vadhan, 2016, §1.6, p. 9].4 The parameter on the other hand should be small but not arbitrarily small. We may think of ≤ 0.01, ≤ 0.1 or ≤ 1 [Dwork et al., 2011, §1], [Dwork and Roth, 2014, §3.5.2, p. 52]. An important property of differential privacy is that it composes Dwork and Roth [2014]. Theorem 1 (Basic composition). If M1 , . . . , Mk are each (, δ)-differentially private then M = (M1 , . . . , Mk ) is (k, kδ)-differentially private. The above is sometimes referred to as (basic) sequential composition as opposed to parallel composition defined next. Theorem 2 (Parallel composition McSherry [2009]). Let Mi each provide (, δ)-differential privacy. Let Xi be arbitrary disjoint subsets of the domain X . The sequence of Mi (D ∩ Xi ) provides (, δ)-differential privacy, where D ∈ X n . A more advanced √ form of composition, which we shall use in this paper, only depletes the “privacy budget” by a factor of ≈ k. See [Dwork and Roth, 2014, §3.5.2] for a precise definition of adaptive composition. Theorem 3 ((Advanced) adaptive composition Dwork et al. [2010], Dwork and Roth [2014], Gaboardi et al. [2014]). Let M1 , . . . , Mk be a sequence of mechanisms where each can take as input the output of a previous mechanism, and let each be 0 -differentially private. Then M(D) = (M1 (D), . . . , Mk (D) is -differentially private for = k0 , and (, δ)-differentially private for r 1 0 = 2k ln 0 + k0 (e − 1), δ for any δ ∈ (0, 1). Furthermore, if the mechanisms Mi are each (0 , δ 0 ) differentially private, then M is (, kδ 0 + δ) differentially private with the same as above. Another key feature of differential privacy is that it is immune to post-processing. Theorem 4 (Post-processing Dwork and Roth [2014]). If M : X n → R is (, δ)-differentially private and f : R → R0 is any randomized function, then f ◦ M : X n → R0 is (, δ)-differentially private. A query is defined as a function q : X n → R. 3 For instance, in D , we may have two binary attributes set to 1, which correspond to two different values of the same B attribute in the original dataset D, e.g., USA and France. 4 A function f in n is negligible, if for all c ∈ N, there exists an n ∈ N such that for all n ≥ n , it holds that f (n) < n−c . 0 0 4

Definition 2 (Counting queries and point functions). A counting query is specified by a predicate q : X → {0, 1} and extended to datasets D ∈ X n by summing up the predicate on all n rows of the dataset as X q(D) = q(x). x∈D A point function Vadhan [2016] is the sum of the predicate qy : X → {0, 1}, which evaluates to 1 if the row is equal to the point y ∈ X and 0 otherwise, over the dataset D. Note that computing all point functions, i.e., answering the query qy (D) for all y ∈ X , amounts to computing the histogram of the dataset D. Definition 3 (Global sensitivity Dwork and Roth [2014], Vadhan [2016]). The global sensitivity of a counting query q : X n → N is ∆q = max 0 n |q(D) − q(D0 )|. D,D ∈X D∼D 0 Definition 4 ((α, β)-utility). The mechanism M is said to be (α, β)-useful for the query class Q if for any q ∈ Q, P [util(M(q, D), q(D)) ≤ α] ≥ 1 − β, where the probability is over the coin tosses of M, and util is a given metric for utility. The Laplace mechanism is employed as a building block in our algorithm to generate the synthetic dataset. Definition 5 (Laplace mechanism Dwork et al. [2006]). The Laplace distribution with mean 0 and scale b has the probability density function 1 |x| Lap(x | b) = e− b . 2b We shall remove the argument x, and simply denote the above by Lap(b). Let q : X n → R be a query. The mechanism MLap (q, D, ) = q(D) + Y where Y is drawn from the distribution Lap (∆q/) is known as the Laplace mechanism. The Laplace mechanism is -differentially private [Dwork and Roth, 2014, §3.3]. Furthermore, with probability at least 1 − β [Dwork and Roth, 2014, §3.3] ∆q 1 . max |q(D) − MLap (q, D, )| ≤ ln = α, q∈Q β where β ∈ (0, 1]. If q is a counting query, then ∆q = 1, and the error α is 1 1 α = ln β with probability at least 1 − β. 2.4 Overview of Copulas For illustration, we assume the bivariate case, i.e., we have only two variables (attributes) X1 and X2 . Let F1 and F2 be their margins, i.e., cumulative distribution functions (CDFs), and let H be their joint distribution function. Then, according to Sklar’s theorem Sklar [1959], Nelsen [2006], there exists a copula C such that for all X1 , X2 ∈ R, H(X1 , X2 ) = C(F1 (X1 ), F2 (X2 )), which is unique if F1 and F2 are continuous; otherwise it is uniquely determined on Ran(F1 ) × Ran(F2 ), where Ran(·) denotes range. In other words, there is a function that maps the joint distribution function to 5

each pair of values of its margins. This function is called a copula. Since our treatment is on binary variables X ∈ {0, 1}, the corresponding margins are defined over X ∈ R as 0, X ∈ [−∞, 0) F (X) = a, X ∈ [0, 1) 1, X ∈ [1, ∞] where, a ∈ I. The above satisfies the definition of a distribution function [Nelsen, 2006, §2.3]. This allows us to define the quasi-inverse of F , denoted F −1 , (tailored to our case) as ( 0, t ∈ [0, a] F −1 (t) = 1, t ∈ (a, 1] Now, using the quasi-inverses, we see that C(u, v) = H(F1−1 (u), F2−1 (v)), where u, v ∈ I. If an analytical form of the joint distribution function is known, the copula can be constructed from the expressions of F1−1 and F2−1 (provided they exist). This then allows us to generate random samples of X1 and X2 by first sampling a uniform u ∈ I, and then extracting v from the conditional distribution. Using the inverses we can extract the pair X1 , X2 . See [Nelsen, 2006, §2.9] for more details. However, in our case, and in general for real-world datasets, we seldom have an analytical expression for H, which could allow us to construct C. There are two approaches to circumvent this, which rely on obtaining an empirical estimate of H from DB . • The first approach relies on constructing a discrete equivalent of an empirical copula [Nelsen, 2006, §5.6, p. 219], using the empirical H obtained from the dataset DB . However, doing this in a differentially private manner requires computing answers to the set of all point functions (Definition 2) of the original dataset Vadhan [2016], which amounts to finding the histogram of DB . Unfortunately, for high dimensional datasets, existing efficient approaches of differentially private histogram release Vadhan [2016] would release a private dataset that discards most rows of the original dataset. This is due to the fact that high dimensional datasets are expected to have high number of rows with low multiplicity. • The other approach, and indeed the one taken by us, is using some existing copula and adjusting its parameters according to the dataset DB . For this paper we choose the Gaussian copula, i.e., the copula C(u, v) = Φr (Φ−1 (u), Φ−1 (v)), where Φ is the standard univariate normal distribution, r is the Pearson correlation coefficient and Φr is the standard bivariate normal distribution with correlation coefficient r. We can then replace the Φ’s with the given marginals F and G resulting in a distribution which is not standard bivariate, if F and G are not standard normal themselves. The underlying assumption in this approach is that the given distribution H is completely characterised by the margins and the correlation. This in general may not be true of all distributions H. However, in practice, this can potentially provide good estimates of one way margins and a subset of the two way margins as discussed earlier. Another advantage of this approach is computational efficiency: the run time is polynomial in m (the number of attributes) and n (the number or rows). While our introduction to copulas has focused on two variables, this can be extended to multiple vari- ables Sklar [1959]. 3 Proposed Mechanism Our method has three main steps as shown in Figure 1: data pre-processing, differentially private statistics, and copula construction. Among these, only the second step involves privacy treatment. The last step preserves privacy due to the post-processing property of differential privacy (Theorem 4). The first step, if not undertaken carefully, may result in privacy violations, as we shall discuss shortly. We shall elaborate each step in the following. 6

Start Data Pre-processing Differentially Private Statistics Copula Construction Obtaining Gaussian Privately Computing Attribute Binning Correlation from Synthetic Dataset One-way Margins Correlation Matrix Generating Rows Binary Privately Computing End through Gaussian Representation Correlations Copula Figure 1: Flowchart of our method. 3.1 Data Pre-processing 3.1.1 Binning of Continuous and Ordinal Data Since we Pwill convert the original data into its binary format we need to ensure that the resulting expansion, m i.e., d = i=1 |Ai |, does not become prohibitively large. This will in general be the case with continuous attributes or discrete valued attributes with large domains. To overcome this, we bin these attributes into discrete intervals. However, care must be taken to ensure that the resulting binning is done independent of the input dataset D to avoid privacy violations. For instance, data-specific binning could introduce a bin due to an outlier which would reveal the presence and absence of that particular individual, violating differential privacy. We, therefore, assume that this binning is done only through what is publicly known about each attribute. For instance, when binning the year of birth into bins of 10 years, we assume that most individuals in the dataset have an age less than 120 (without checking this in the dataset). This of course depends on the nature of the attributes. While binning can be done in a differentially private manner, we use a simpler approach as this is not the main theme of our paper. We also note that data binning as a pre-processing step is almost always adopted in related work, e.g., Gaboardi et al. [2014]. 3.1.2 Binary Representation In order to use the Gaussian copula, we need to determine correlations between attributes via some fixed correlation coefficient; which is a numeric measure. We use the Pearson product moment correlation co- efficient for this purpose. This means that in order to measure correlations, either the attribute values in the dataset need to be real valued or have to mapped to a real number. These mapped values then follow their order in the set of real numbers. If an attribute is ordinal, e.g., age or salary, a natural order exists. However, categorical attributes do not have a natural order, e.g., gender. One way to solve this conundrum is to induce an artificial order among the values of the categorical attribute, e.g., through a hierarchical tree Iyengar [2002]. However, this ordering is inherently artificial and any change in order results in a markedly different correlation. An illustration of this fact is shown in Appendix A (we also briefly discuss this in Section 5).5 Our approach instead is to work on the binary version of the dataset obtained via dummy coding, as explained in Section 2.2. Pairwise correlations now amount to finding Pearson correlation between pairs of binary variables. This way of interpreting data makes no distinction between ordinal and categorical variables by not assigning any order to the latter. 5 There are other correlation measures that can be used for categorical attributes, e.g., Cramér’s V Cramér [2016]. However, we need a unified measure that is applicable to all attributes alike. 7

3.2 Differentially Private Statistics 3.2.1 Privately Computing One-Way Margins The first ingredient in our method is the one-way margins, i.e., marginal CDFs, of the attributes Xi , for i ∈ {1, 2, . . . , d}. We denote these margins by F̂i (x), where x ∈ {0, 1}. Since each Xi is binary, these margins can be calculated as Pn (k) I{Xi = 0} F̂i (0) = k=1 , F̂i (1) = 1 n Pn (k) To make F̂i (x) differentially private, we add Laplace noise to the counts n̂0i = k=1 I{Xi = 0} and n (k) n̂1i k=1 I{Xi = 1} to obtain F̃i (x), which is summarized in Algorithm 1. If the differentially private sum P is negative, we fix it to 0. Note that this utilizes the post-processing property (cf. Theorem 4) and hence maintains privacy. The reason we add noise to both n̂0i and n̂1i instead of just adding noise to n̂0i is to avoid the noisy n̂0i exceeding n. ALGORITHM 1: Obtaining Differentially Private Marginals F̃i (x) 1. For the ith attribute, count the numbers of events Xi = 0 and Xi = 1 as n X n o n X n o (k) (k) n̂0i = I Xi = 0 , n̂1i = I Xi = 1 k=1 k=1 2. Add Laplace noise to n̂0i and n̂1i , and obtain the noisy counts as 2 2 ñ0i = n̂0i + Lap 0 , ñ1i = n̂1i + Lap 0 i i where 0i is the privacy budget associated with computing the ith margin. 3. Obtain F̃i (x) as ñ0i F̃i (0) = , F̃i (1) = 1 ñ0i + ñ1i Algorithm 1 is (0i , 0)-differentially private due to differential privacy of the Laplace mechanism and the fact that the algorithm is essentially computing the histogram associated with attribute Xi [Dwork and Roth, 2014, §3.3, p. 33]. Importantly, the privacy budget 0i is impacted only by the number of attributes m in the original dataset D and not by the number of binary attributes d in the binary version DB . This is shown by the following lemma. Lemma 1. Fix an attribute A in D, and let X1 , . . . , X|A| denote the binary attributes constructed from A. Then if the computation of each marginal Xj , j ∈ [|A|], is (0i , 0)-differentially private, the computation of the marginal A is (0i , 0)-differentially private. Proof. See Appendix B. 3.2.2 Privately Computing Correlations The other requirement of our method is the computation of pairwise correlations given by E {Xi Xj } − µi µj r̂i,j = p , i, j ∈ {1, 2, . . . , d}. (1) var (Xi ) var (Xj ) To obtain the differentially private version of r̂i,j , denoted r̃i,j , one way is to compute r̂i,j directly from DB and then add Laplace noise scaled to the sensitivity of r̂. However, as we show in Appendix C, the empirical correlation coefficient from binary attributes has high global sensitivity, which would result in highly noisy r̃i,j . The other approach is to compute each of the terms in Eq. 1 in a differentially private manner and 8

then obtain r̃i,j . Notice that for binary attribute Xi , its mean µi is given by n̂1i /n. This can be obtained differentially privately as ñ0 µ̃i = 1 − F̃i (0) = 1 − 0 i 1 , ñi + ñi which we have already computed. Likewise, the variance var(X ˆ i ) is given by µ̂i (1 − µ̂i ), whose differentially private analogue, i.e., var(X ˜ i ), can again be obtained from the above. Thus, the only new computation is Pn (k) the computation of E {Xi Xj }, which for binary attributes is equivalent to computing n1 k=1 I{(Xi = (k) 1) ∧ (Xj = 1)}. Algorithm 2 computes this privately. ALGORITHM 2: Obtaining Differentially Private Two-Way Positive Conjunctions Ẽ {Xi Xj } 1. For the ith and jth attributes (1 ≤ i < j ≤ d), and for a, b ∈ {0, 1}, count the number of events (Xi , Xj ) = (a, b) as n X n o (k) (k) n̂ab i,j = I (Xi = a) ∩ (Xj = b) (2) k=1 2. Add Laplace noise onto n̂ab i,j for a, b ∈ {0, 1} and obtain the noisy count as 2 ñab ab i,j = n̂i,j + Lap , (3) 00i,j where 00i,j is the privacy budget associated with computing the (i, j)th two-way marginal. 3. Obtain Ẽ {Xi Xj } as ñab i,j Ẽ {Xi Xj } = P . (4) a,b∈{0,1} ñab i,j Algorithm 2 is (00i,j , 0)-differentially private due to the differential privacy of the Laplace mechanism and the fact that the algorithm is essentially computing the histogram associated with attribute pairs (Xi , Xj ) [Dwork and Roth, 2014, §3.3, p. 33]. Once again, the privacy budget 00i,j is impacted only by the number of pairs of attributes m 2 in the original dataset D and not by the number of pairs of binary attributes d2 in the binary version DB . This is presented in the following lemma, whose proof is similar to that of Lemma 1, and hence omitted for brevity. Lemma 2. Fix two attributes Ai and Aj , i 6= j, in D. Let Xi,1 , . . . , Xi,|Ai | and Xj,1 , . . . , Xj,|Aj | denote the binary attributes constructed from Ai and Aj , respectively. Then, if the computation of each of the two-way marginals (Xi,k , Xj,k0 ), k ∈ [|Ai |], k 0 ∈ [|Aj |], is (00i,j , 0)-differentially private, the computation of all all two-way marginals of Ai and Aj is (00i,j , 0)-differentially private. The differentially private correlation coefficients r̃i,j thus obtained can be readily used to construct the differentially private correlation matrix r̃1,1 r̃1,2 ··· r̃1,d r̃2,1 r̃2,2 ··· r̃2,d R̃ = . . (5) .. .. .. .. . . . r̃d,1 r̃d,2 ··· r̃d,d Notice that the above algorithm computes the correlations in time O(d2 n). This can be prohibitive if d is large, i.e., if each of the m (original) attributes have large domains. An alternative algorithm to compute the two-way positive conjunctions that takes time only O(m2 n) is shown in Appendix D. 3.3 Copula Construction For this section, we do not need access to the database any more. Hence, any processing done preserves the previous privacy budget due to closure under post-processing (Theorem 4). 9

3.3.1 Obtaining Gaussian Correlation from the Correlation Matrix Our aim is to sample standard normal Gaussian variables Yi ’s corresponding to the attributes X̃i ’s6 where the correlations among Yi ’s, given by the correlation matrix P, are mapped to the correlations among the Xi ’s, given by the (already computed) correlation matrix R̃. A sample from the Gaussian variable Yi is transformed backed to X̃i as X̃i = F̃i−1 (Φ (Yi )) . (6) Obviously, if the attributes X̃i are independent, then Yi are also independent. However, in practice X̃i ’s are correlated, which is characterized by R̃ in our case. Hence, the question becomes: How to choose a correlation matrix P for Yi ’s, so that X̃i ’s have the target correlation relationship defined by R̃? From Eq. 1 and its perturbation through r̃i,j , we can obtain n o r E X̃i X̃j = r̃i,j var X̃i var X̃j + µ̃i µ̃j , (7) On the other hand, from Eq. 6, we can get n o n o E X̃i X̃j = E F̃i−1 (Φ (Yi )) F̃j−1 (Φ (Yj )) Z +∞ Z +∞ = F̃i−1 (Φ (yi )) F̃j−1 (Φ (yj )) Φi,j (yi , yj ) dyi dyj , (8) −∞ −∞ where Φi,j (yi , yj ) denotes the standard bivariate probability density function (PDF) of the correlated stan- dard normal random variables Yi and Yj , given by ! 1 yi2 + yj2 − 2ρi,j yi yj Φi,j (yi , yj ) = exp − . (9) 2 1 − ρ2i,j q 2π 1 − ρ2 i,j Here ρi,j is the Pearson correlation coefficient, which for the standard normal variables Yi and Yj is given by E {Yi Yj }. Our task is to find the value of ρi,j such that Eq. 7 and Eq. 8 are equal. In other words, Eqs. 7 and 8 define the relationship between r̃i,j and ρi,j . Notice that by construction, in general, we do not have ρi,j = r̃i,j . We can obtain ρi,j by means of a standard bisection search (e.g., see [Nelson, 2015, p. 148]). The two-fold integral in Eq. 8 with respect to yi and yj , is evaluated numerically in the bisection search. In more detail, dyi and dyj are set to a small value of 0.01, and the lower and upper limits of the such integral are set to -10 and 10, respectively. Such lower and upper limits make the numerical results sufficiently accurate since Φi,j (yi , yj ) is the standard bivariate PDF of two correlated normal random variables (and hence have negligibly small probability mass beyond the limits of integration). With each such ρi,j , we can construct the matrix P corresponding to R̃. Remark. Notice that the choice of ρi,j ensures that the resulting n sampled o Gaussian variables have the property that when transformed back to X̃i and X̃j , we get E X̃i X̃j ≈ E {Xi Xj }, where the latter is n o Pn (k) (k) the input expectation. For binary attributes, recall that E {Xi Xj } = n1 k=1 I (Xi = 1) ∩ (Xj = 1) . Thus the method ensures that the “11”’s are well approximated to the input distribution. However, the correlation coefficient, and the distribution of 01’s, 10’s and 00’s might not be the same as the original dataset. This is evident from Eq. 7. However, since the remaining quantities in Eq. 7 depend on the one-way margins, maintaining a good approximation to margins would imply that the distribution of 01’s, 00’s and 10’s would also be well approximated. However, error in one or both of the margins would propagate to the error in the 01’s, 10’s and 00’s. It is due to this reason that we target positive two-way conjunctions for utility. 6 i.e., the synthetic versions of the Xi ’s. 10

ALGORITHM 3: Algorithm to Obtain NCM (P) 1. Initialization: 4S0 ← 0,Y0 ← P, k ← 1. 2. While k < iters, where iters is the maximum of iterations (a) Rk = Yk−1 − 4Sk−1 (b) Projection of Rk to a semi-positive definite matrix: Xk ← VT diag (max {Λ, 0}) V, where V and Λ contain the eigenvectors and the eigenvalues of Rk , respectively, and diag (·) transforms a vector into a diagonal matrix. (c) 4Sk ← Xk − Rk (d) Projection of Xk to a unit-diagonal matrix: Yk ← unitdiag (Xk ), where unitdiag(·) fixes the diagonal of the input matrix to ones. (e) k ← k + 1. 3. Output: P0 ← Yk . 3.3.2 Generating Records Following a Multivariate Normal Distribution Since each Yi follows the standard normal distribution and the correlation among Yi ’s is characterized by P, we can generate records by sampling from the resulting multivariate normal distribution. However, the matrix P obtained through this process should have two important properties for this method to work: 1. P should be a correlation matrix, i.e., P should be a symmetric positive semidefinite matrix with unit diagonals Higham [2002]. 2. P should be positive definite to have a unique Cholesky decomposition defined as P = LT L where L is a lower triangular matrix with positive diagonal entries [Golub and Van Loan, 2013, p. 187]. To ensure Property 1, we use the algorithm from [Higham, 2002, §3.2], denoted NCM, to obtain the matrix P0 as P0 = NCM (P) , (10) Furthermore, to satisfy Property 2, i.e., positive definiteness, we force the zero eigenvalues of P0 to small positive values. For completeness, we describe the simplified skeleton of the algorithm from Higham [2002] in Algorithm 3. The algorithm searches for a correlation matrix that is closest to P in a weighted Frobenius norm. The output is asymptotically guaranteed to output the nearest correlation matrix to the input matrix Higham [2002]. For details see [Higham, 2002, §3.2, p. 11]. The overall complexity of the procedure is O(dω ), where ω < 2.38 is the coefficient in the cost of multiplying two d × d matrices Bardet et al. [2003]. Note that matrix multiplication is highly optimized in modern day computing languages. Finally, we generate records from a multivariate normal distribution using the well known method de- scribed in Algorithm 4. The rationale of invoking Cholesky decomposition is to ensure that E YT Y = E LT ZT ZL = LT E ZT Z L = LT L = P0 , where we have used the fact that E ZT Z = I because each record in Z follows i.i.d. multivariate normal distribution. The output dataset DB0 is the final synthetic dataset. 3.4 Privacy of the Scheme Let m = m + m 0 2 , where m is the number of attributes in D. Let < 1, say = 0.99. Fix a δ. We set each of the 0i ’s and 00i,j ’s for i, j ∈ [m] to 0 , where 0 is such that it gives us the target through Theorem 3 by setting k = m0 . According to the advanced composition theorem (Theorem 3), since each of our m0 mechanisms are (0 , 0) differentially private, the overall construction is (, δ)-differentially private. Privacy budget consumed over each of the m0 mechanisms is roughly √1m0 . Note that all the algorithms in Section 3.3 do not require access to the original dataset and therefore privacy is not impacted due to the post-processing property (see Theorem 4). 11

ALGORITHM 4: Generating Records Following a Multivariate Normal Distribution 1. Generate n records, with the attributes in each record following i.i.d. standard normal distribution, i.e., T Z = Z (1) Z (2) Z (n) ··· h iT (k) where Z (k) is given by Z1(k) Z2 (k) ··· Zd (k) and Zi ’s are i.i.d. standard normal random variables. 2. Compute Y as Y = ZL where LT L = P0 and L is obtained by the Cholesky decomposition as L = chol (P0 ). T (k) (k) 3. From the matrix Y = Y (1) Y (2) · · · Y (n) , map every element Yi in each record Y (k) to X̃i using Eq. 6, where i ∈ {1, 2, . . . , d}. 4. Output the mapped data as DB0 . 4 Experimental Evaluation and Utility Analysis 4.1 Query Class and Error Metrics As mentioned in Section 1 and explained in Section 3.3.1, our focus is on all one-way marginal and two-way positive conjunction queries. In addition, we will also evaluate the performance of our method on the set of three-way positive conjunction queries to demonstrate that our method can also give well approximated answers to other types of queries. Thus, we use the following query class to evaluate our method, Q := Q1 ∪ Q2 ∪ Q3 . Here, Q1 is the set of one-way marginal counting queries, which consists of queries q specified as n (k) X q(i, DB ) = I{Xi = b} k=1 where i ∈ [d] and b ∈ {0, 1}. The class Q2 is the set of positive two-way conjunctions and consists of queries q specified as n (k) (k) X q(i, j, DB ) = I{(Xi = 1) ∩ (Xj = 1)} k=1 where i, j ∈ [d], j 6= i. We define Q12 = Q1 ∪ Q2 . The class Q3 is the set of positive three-way conjunctions, and is defined analogously. Note that only those queries are included in Q2 and Q3 whose corresponding binary columns are from distinct original columns in D. Answers to queries which evaluate at least two binary columns from the same original column in D can be trivially fixed to zero; as these are “structural zeroes.” We assume this to be true for queries in the two aforementioned query sets from here onwards. Our error metric of choice is the absolute error, which for a query q ∈ Q is defined as |q(DB ) − q(DB0 )| We have preferred this error metric over relative error (which scales the answers based on the original answer), since it makes it easier to compare results across different datasets and query types. For instance, in the case of relative error, a scaling factor is normally introduced, which is either a constant or a percentage of the number of rows in the dataset Xiao et al. [2011], Li et al. [2014]. The scaling factor is employed to not penalize queries with extremely small answers. However, the instantiation of the scaling factor is mostly a heuristic choice. For each of the datasets (described next), we shall evaluate the differences in answers to queries from Q in the remainder of this section. We shall be reporting (α, β)-utility in the following way. We first sort the query answers in ascending order of error. We then fix a value of β, and report the maximum error and average error from the first 1 − β fraction of queries. We shall use values of β = 0.05, 0.01 and 0. The errors returned then correspond to 95%, 99% and 100% (overall error) of the queries, respectively. 12

Differential Output Mechanism Synthetic Privacy cop Gaussian copula with original correlation matrix 7 3 cop-ID Gaussian copula with identity correlation matrix 7 3 cop-1 Gaussian copula with correlation matrix of all 1’s 7 3 no-cor Answers generated assuming no correlation 7 7 Lap Laplace mechanism 3 7 dpc Gaussian copula with diff. private correlation matrix 3 3 Table 1: Notation used for outputs from different mechanisms. Note that not all of them are synthetic datasets or differentially private. 4.2 Datasets and Parameters We used three real-world datasets to evaluate our method. All three datasets contained a mixture of ordinal (both discrete and continuous) and categorical attributes. 1. Adult Dataset: This is a publicly available dataset which is an extract from the 1994 US census information Lichman [2013]. There are 14 attributes in total which after pre-processing result in 194 binary attributes. There a total of 32,560 rows in this dataset (each belonging to an individual). 2. DSS Dataset: This dataset is a subset of a dataset obtained from the Department of Social Services (DSS), a department of the Government of Australia. The dataset contains transactions of social security payments. The subset contains 27 attributes, which result in 674 binary attributes. There are 5,240,260 rows in this dataset. 3. Hospital Dataset: This dataset is a subset of the hospital ratings dataset extracted from a national patient survey in the US.7 The dataset is among many other datasets gathered by the Centers for Medicare & Medicaid Services (CMS), a federal agency in the US. We extracted (the first) 9 attributes and 10,000 rows for this dataset (resulting in 1,201 binary attributes). We mainly chose this dataset as an example of a dataset with highly correlated attributes. Parameter Values. We set δ = 2−30 , following a similar value used in [Gaboardi et al., 2016, §3, p. 5]. This is well below n−1 for the three datasets, where n is the number of rows. We set the same privacy budget, 0 , to compute each of the m one-way marginals and m 2 two-way positive conjunctions. For the Adult, DSS and Hospital datasets we have m = 14, m = 27 and m = 9, respectively. We search for an 0 for the two datasets by setting δ = 2−30 and k = m + m 2 in Theorem 3 which gives an overall of just under 1. The resulting computation gives us 0 = 0.014782 for Adult, 0 = 0.007791 for DSS, and 0 = 0.022579 for the Hospital dataset. 4.3 Experimental Analysis To evaluate our method, we generate multiple synthetic datasets from each of the three datasets. We will first evaluate the synthetic dataset generated through the Gaussian copula (with no differential privacy) for each of the three datasets. This will be followed by the evaluation of the differentially private version of this synthetic dataset against the baseline which is the application of the Laplace mechanism on the set of queries Q. Note that this does not result in a synthetic dataset. For readability, we use abbreviations for the different outputs. These are shown in Table 1. 4.3.1 Error due to Gaussian Copula without Differential Privacy We first isolate and quantify query errors from a synthetic dataset obtained directly through the Gaussian copula, i.e., without any differentially private noise added to the one-way marginals and the correlation matrix (see Figure 1). Since generating synthetic datasets through the copula is an inherently random 7 See https://data.medicare.gov/Hospital-Compare/Patient-survey-HCAHPS-Hospital/dgck-syfz. 13

process, this itself may be a source of error. We denote such a dataset by “cop.” Thus, Adult cop, DSS cop and Hospital cop are the “cop” versions of the corresponding datasets. We restrict ourselves to the query set Q12 and compare error from cop against three other outputs: cop-ID: A synthetic dataset obtained by replacing the correlation matrix P0 with the identity matrix. Eval- uating against this dataset will show whether cop performs better than a trivial mechanism which assumes all binary attributes to be uncorrelated. Note that this mainly effects answers to Q2 , and not the one-way marginals Q1 . cop-1: Another synthetic dataset obtained by replacing the correlation matrix P0 with the matrix 1 of all ones. This serves as other extreme where all attributes are assumed to be positively correlated. Once again, this is to compare answers from Q2 . no-cor: For the query class Q2 , we obtain a set of random answers which are computed by simply multiplying the means µ1 and µ2 of two binary attributes. This is the same as simulating two-way positive conjunctions of uncorrelated attributes. This should have an error distribution similar to the answers on Q2 obtained from cop-ID. Results. Figure 2 shows the CDF of the absolute error on the query set Q12 from different outputs from the three datasets. Looking first at the results on Q1 (top row in the figure), we see that for all three datasets cop has low absolute error, yielding (α, β)-utility of (149, 0.01) for Adult, (170, 0.01) for DSS and (28, 0.01) for the Hospital dataset in terms of max-absolute error. The datasets cop-ID and cop-1 exhibit similar utility for all three datasets. This is not surprising since Q1 contains one-way marginals, whose accuracy is more impacted by the inverse transforms (Eq. 6) rather than the correlation matrix. Answers on Q2 are more intriguing (bottom row of Figure 2). First note that the utility from cop is once again good for all three input datasets with 99% of the queries having a maximum absolute error of 353 for Adult, 83 for DSS and 6 for the Hospital dataset. The cop-1 outputs have poorer utility. However, interestingly, cop-ID and no-cor outputs yield utility very similar to cop. We discuss this in more detail next. 1.0 1.0 1.0 0.8 0.8 0.8 Cumulative Distribution Cumulative Distribution Cumulative Distribution 0.6 0.6 0.6 0.4 0.4 0.4 Adult cop Dss cop 0.2 Hospital cop 0.2 0.2 Adult cop-ID Dss cop-ID Hospital cop-ID Adult cop-1 Dss cop-1 0.0 Hospital cop-1 0 50 100 150 200 250 0 50 100 150 200 250 0 20 40 60 80 100 120 Absolute Error Absolute Error Absolute Error (a) Set Q1 (one-way marginals) (b) Set Q1 (one-way marginals) (c) Set Q1 (one-way marginals) 1.0 1.0 1.0 0.8 0.8 0.8 Cumulative Distribution Cumulative Distribution Cumulative Distribution 0.6 0.6 0.6 0.4 0.4 0.4 0.2 adult cop 0.2 Dss cop Hospital cop adult cop-ID Dss cop-ID 0.2 Hospital cop-ID adult cop-1 Dss cop-1 Hospital cop-1 0.0 adult no-cor 0.0 Dss no-cor 0.0 Hospital no-cor 0 200 400 600 800 1000 0 200 400 600 800 1000 0 100 200 300 400 500 Absolute Error Absolute Error Absolute Error (d) Set Q2 (two-way conjunctions) (e) Set Q2 (two-way conjunctions) (f) Set Q2 (two-way conjunctions) Figure 2: Relative error over Q12 on synthetic versions of the Adult, DSS and Hospital datasets (without differential privacy). 14

Separating Low and High Correlations. There are two possible reasons why cop-ID and no-cor outputs perform close to cop outputs on the set Q2 : (a) our method does not perform better than random (when it comes to Q2 ), or (b) uncorrelated attributes dominate the dataset thus overwhelming the distribution of errors. The second reason is also evidenced by the fact that the cop-1 dataset, using a correlation matrix of all ones, performs worse than the other three outputs, indicating that highly correlated attributes are rare in the three datasets. To further ascertain which of the two reasons is true, we separate the error results on the set Q2 into two parts: a set of (binary) attribute pairs having high correlations (positive or negative), and another with low correlations. If our method performs better than random, we would expect the error from cop to be lower on the first set when compared to cop-ID and no-cor, while at least comparable to the two on the second set. We use the Hospital dataset for this analysis as it was the dataset with the most highly correlated attributes. There are a total of 626,491 pairs of binary attributes (as mentioned before, we ignore binary attribute pairs that correspond to different attribute values on the same attribute in the original dataset). Out of these, only 3,355 have an (absolute) Pearson correlation coefficient |r| ≥ 0.5. Thus, an overwhelming 99.46% of pairs have |r| < 0.5. This shows that the dataset does indeed contain a high number of low-correlated attributes, which partially explains similar error profile of cop, cop-ID and no-cor. Figure 3 shows this breakdown. The error on the set with |r| < 0.5 is very similar for cop, cop-ID and no-cor (Figure 3(a)). Looking at Figure 3(b), for the set with |r| ≥ 0.5, on the other hand, we note that cop outperforms both cop-ID and no-cor. Also note that cop-1 is similar in performance to our method. This is understandable since cop-1 uses a correlation matrix with all ones, and hence is expected to perform well on highly correlated pairs. This indicates that our method outperforms cop-ID and no-cor. We conclude that the apparent similarity between cop, cop-ID and no-cor on the set Q2 is due to an artefact of some real-world datasets which may have a high number of uncorrelated attributes; when the results are analyzed separately, our method is superior. We remark that we arrived at the same conclusion for the Adult and DSS datasets, but omit the results due to repetition. Pearson R < 0.5 (623136/626491) 3355/626491 pairs of attributes 1.0 1.0 0.8 0.8 Cumulative Distribution Cumulative Distribution 0.6 0.6 0.4 0.4 Hospital cop 0.2 Hospital cop 0.2 Hospital cop-ID Hospital cop-ID Hospital cop-1 Hospital cop-1 0.0 Hospital no-cor 0.0 Hospital no-cor 0 50 100 150 200 250 0 50 100 150 200 250 Absolute Error Absolute Error (a) |r| < 0.5 (b) |r| ≥ 0.5 Figure 3: Absolute error over Q2 on synthetic versions of the Hospital dataset (without differential privacy) for different values of the absolute Pearson correlation coefficient |r| between pairs of attributes. 4.3.2 Error due to Differentially Private Gaussian Copula Having established that the error due to Gaussian copula is small, we now turn to the complete version of our method, i.e., with differential privacy. For this section, we are interested in three outputs: (a) dpc, i.e., the synthetic dataset obtained through our differentially private Gaussian copula method, (b) cop, i.e., the synthetic dataset via our Gaussian copula method without differential privacy, and (c) Lap, i.e., a set of answers obtained by adding independent Laplace noise to the answers to the queries in Q12 . Note that Lap is not a synthetic dataset. We set the same 0 for Lap, as we did for dpc. Results. Figure 4 shows the absolute error CDF on the query set Q12 for cop, dpc and Lap versions constructed from the three datasets. For the set Q1 (top row in figure), we see that cop outperforms both dpc and Lap. This is due to the fact that for privacy, a higher amount of noise is required. Crucially, 15

our method does not introduce further error over the Laplace mechanism, as is indicated by the similarity of the curves corresponding to dpc and Lap. Interestingly, the results for Q2 show that dpc outperforms independent Laplace noise (bottom row of Figure 4). While the error from dpc is still higher than cop, it is closer to it than the error due to Lap. This indicates that for the majority of the queries, our method applies less noise than Lap. However, in some cases, for a small percentage of queries Lap adds less noise than our mechanism. This is clear from Table 2, where we show the maximum and average error8 from 95%, 99% and 100% percent of the queries from Q12 across dpc and Lap versions of the Adult, DSS and Hospital datasets. The error profiles of the dpc and Lap variants are similar for the query set Q1 . For the set Q2 , we can see that Lap only outperforms our method for the Adult and DSS datasets if we consider the maximum absolute error across all queries. On the other hand our method outperforms Lap if we consider 95% and 99% of queries. Thus, for less than 1% of queries, the dpc versions of Adult and DSS exhibits less utility than Lap, while being similar to the latter for the Hospital dataset. 1.0 1.0 1.0 0.8 0.8 0.8 Cumulative Distribution Cumulative Distribution Cumulative Distribution 0.6 0.6 0.6 0.4 0.4 0.4 0.2 Adult cop 0.2 Dss cop 0.2 Hospital cop Adult dpc Dss dpc Hospital dpc 0.0 Adult Lap 0.0 Dss Lap 0.0 Hospital Lap 0 100 200 300 400 500 600 700 0 200 400 600 800 1000 1200 1400 0 100 200 300 400 500 Absolute Error Absolute Error Absolute Error (a) Set Q1 (Adult) (b) Set Q1 (DSS) (c) Set Q1 (Hopspital) 1.0 1.0 1.0 0.9 0.9 0.9 Cumulative Distribution Cumulative Distribution Cumulative Distribution 0.8 0.8 0.7 0.8 0.7 0.6 0.7 0.6 0.5 0.5 0.6 0.4 Adult cop 0.4 Dss cop Hospital cop 0.3 Adult dpc Dss dpc 0.5 Hospital dpc 0.2 Adult Lap 0.3 Dss Lap Hospital Lap 0 200 400 600 800 1000 1200 1400 0 200 400 600 800 1000 1200 1400 0 200 400 600 800 1000 Abolute Error Abolute Error Abolute Error (d) Set Q2 (Adult) (e) Set Q2 (DSS) (f) Set Q2 (Hopsital) Figure 4: Absolute error over Q12 of one-way marginals (set Q1 ) and two-way positive conjunctions (set Q2 ) on synthetic versions of the Adult (left), DSS (middle) and Hospital (right) datasets with differential privacy. 4.3.3 Results on Three-Way Conjunctions Even though our method is expected to perform best on the query set Q12 , we show that the method performs well on other types of queries as well. For this, we use the set of three-way positive conjunction queries as an example, i.e., Q3 . We compare the error against the answers returned from (independent) Laplace mechanism by choosing an appropriate value of 0 for each query in Q3 according to the advanced composition theorem (see Section 4.2), such that overall is just under 1 for queries in the set Q3 only (i.e., the Laplace mechanism does not compute answers to Q12 ). The results are shown in Figure 5. Once again our method outperforms the Laplace mechanism for the majority of the queries in Q3 for all three datasets. Looking closely, we see from Table 2, that Lap actually performs better in terms of maximum absolute error for 1% of the queries in the Adult dataset. However, for majority of the queries, > 99%, over 8 Rounded to the nearest integer. 16

Adult DSS Hospital Mechanism 95% 99% 100% 95% 99% 100% 95% 99% 100% ave max ave max ave max ave max ave max ave max ave max ave max ave max dpc (one-way) 92 389 107 482 106 773 151 657 175 982 185 1267 53 208 61 302 65 622 Lap (one-way) 85 353 99 505 105 741 149 621 172 963 182 1287 52 199 61 299 61 887 dpc (two-way) 18 184 29 504 38 4788 7 108 15 426 23 5744 1 5 1 22 1 836 Lap (two-way) 58 317 72 536 78 1246 99 594 125 995 136 2947 31 199 40 333 43 1119 dpc (three-way) 12 120 20 408 28 6148 9 98 16 372 36 9238 1 3 3 55 3 616 Lap (three-way) 102 591 128 1002 138 2508 268 1651 342 2805 726 9562 45 281 57 477 63 1361 Table 2: Absolute error α of the dpc and Lap mechanisms on the three datasets. The columns show average and maximum values of α for 95%, 99% and 100% of the queries (corresponding to β = 0.05, 0.01 and 0.00, respectively). indicates our method significantly outperforms Lap; indicates Lap significantly outperforms our method; significance is defined as an error ratio of approximately 2 or more. method performs better. In fact, our method outperforms Lap in terms of the 95% and 99% error profiles for all three datasets. For the DSS dataset the maximum error over all queries is similar to Lap, whereas for the Hospital dataset we again outperform Lap. 1.0 1.0 1.0 0.8 0.8 0.8 Cumulative Distribution Cumulative Distribution Cumulative Distribution 0.6 0.6 0.6 0.4 0.4 0.4 0.2 0.2 0.2 Adult dpc Dss dpc Hospital dpc 0.0 Adult Lap 0.0 Dss Lap 0.0 Hospital Lap 0 250 500 750 1000 1250 1500 1750 2000 0 500 1000 1500 2000 2500 3000 0 200 400 600 800 1000 Absolute Error Absolute Error Absolute Error (a) Adult (b) DSS (c) DSS Figure 5: CDFs of the absolute error from our method against the Laplace mechanism on the set Q3 of three-way positive conjunctions on the three datasets. 4.3.4 Effect of the Privacy Parameter To show the effect of on utility, we vary it from 0.25 to 5 and report the error on the set Q12 . For this, we only use the Adult dataset as the effect is similar on the other two datasets. Figure 6 shows the CDF of the absolute error against different values of . Notice that this is the overall privacy budget. With = 0.25 we have average and maximum absolute errors of 357 and 3419, respectively, for the set Q1 , and 68 and 6421, respectively, for the set Q2 . With = 5, the average and maximum absolute errora are much lower at 41 and 179, respectively, for Q1 , and 27 and 5882, respectively, for the set Q2 . As expected, the error profiles gradually improve as we move from = 0.25 to = 5. The error profiles are much similar for the set Q2 then Q1 . As discussed in Section 4.3.1, since the majority of attributes are uncorrelated, this implies that our method maintains that aspect by not adding too much noise on the set Q2 . 4.3.5 Computational Time and Parallelism One of the motivations for using the proposed approach is its computational feasibility when the input data is high dimensional. We first note that our method is highly parallelizable. In particular, the computation of one-way and two-way positive conjunctions can be done in parallel. For the one-way marginals, parallel computation is straightforward. For the two-way conjunctions, we take the ith attribute (in the original dataset) and compute its conjunction with all attributes numbered i + 1 to m attributes, for all i ∈ [1, m − 1], assigning a separate process for each i. Obviously, the number of combinations for the first attribute is the highest, and becomes progressively less for latter attributes. Likewise, we also parallelize the computation of d2 − d Gaussian correlations ρi,j which uses a bisection search. While other components of our method 17

You can also read