Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
132,123 characters · 19 sections · 106 citation commands
Differentially Private Estimation via Statistical Depth
Since dwork2006calibrating first introduced the concept of differential privacy (DP), it has become the gold standard notion of privacy in the statistical disclosure limitation literature, both because of the strong privacy guarantees that DP mechanisms provide and the theoretical properties that make the formulation of many DP mechanisms straightforward, which are described in more detail in the next section. However, DP linear regression methods provide a few interesting theoretical issues related to the lack of exogenous bounds on the regression estimates and/or observations in typical use cases. At the same time, since linear regressions are a ubiquitous tool in applied statistics, there are compelling use cases for linear regressions using sensitive data on respondents in small samples. For example, chetty2018opportunity use data from the US Census Bureau and the Internal Revenue Service to estimate socioeconomic mobility of populations within each Census tract based on ordinary least squares regression (OLS) estimates.
This paper explores the use of two notions of statistical depth to formulate approximate DP mechanisms for estimators of location and linear regression coefficients. Specifically, we formulate approximate DP estimators for the Tukey median, which is a multivariate generalization of the median and is defined as the maximizer of halfspace depth, and for the deepest regression, which is defined as the maximizer of regression depth and is an estimator of the median of the dependent variable conditional on the covariates tukey1975mathematics,rousseeuw1999regression. One break we make from the norm in the DP literature is that the mechanisms proposed here do not require bounds on the observations in the dataset, but instead require bounds on the space of feasible estimates. In the case of measures of central tendency like the Tukey median, this alternative is a strict relaxation of the requirement of specifying a bounded set containing the observations, since a bounded set containing the observations also contains all reasonable estimators of central tendency, but the converse is not true in the generic case.
We also provide variants of these methods that satisfy random differential privacy (RDP), which is a relaxation of differential privacy provided by hall2013random. This definition of privacy protects against accurate inferences on individual observations of samples that are sufficiently likely to be drawn from the same population distribution, without attempting to limit inferences on the population distribution itself. The RDP variants of these estimators have the advantage of not requiring bounds on the estimates or the observations. Using bounds on the estimates themselves in our proposed DP mechanisms is also required for our derivations of the proposed RDP mechanisms, since the RDP mechanisms simply call our proposed approximate DP methods after using the dataset itself to define the feasible sets by nonparametric confidence regions for the estimators. Also, unlike the motivating example provided by hall2013random, for many input datasets, all of the proposed RDP mechanisms provided in this paper are invariant to a privacy attacker's prior information set, in the sense that, for any such information set, these mechanisms do not reveal any respondent's data with certainty for a wide class of input datasets. The final mechanism we introduce is an approximate DP Medsweep mechanism, which is an approximation of the deepest regression provided by rousseeuw1998computing that is more computational efficient, particularly in the multivariate setting.
After introducing these estimators, we provide simulations to compare these methods to existing DP regression techniques. One advantageous feature of the DP estimators proposed here is that their statistical performance is typically less dependent on the input bounds for larger sample sizes than existing approaches, which we explore further in these simulations. Specifically, our implementation choices related to bounds on observations and the estimator are intended to err on the side of understating the relative accuracy of the proposed DP regression estimators. The estimators proposed here appear to perform favorably relative to the other DP regression methods considered when either the sample size is larger than approximately 100 and/or $\epsilon$ is sufficiently high. We also use these simulations to estimate the maximum possible rate at which the diameter of the feasible set of estimates can be increased without having any impact on the distribution of the DP estimators. In these simulations, it appears possible to define this feasible set so that its diameter increases at an exponential rate in the sample size. In contrast, the accuracy of many classical DP location and regression estimators depends strongly on the tightness of bounds on the observations. When these bounds are not provided by the use case at hand, one can estimate them using a preliminary DP mechanism; however, one motivation for the RDP estimators proposed here is that this is not always straightforward. For example, chen2016differentially propose using a preliminary DP mechanism that outputs a bounding box of the form $[-c,c]^d$ such that the proportion of datapoints in $[-c,c]^d$ is approximately equal to the user choice parameter $\psi\in (0,1),$ but this approach is not ideal when the observations are centered around a point that is far from the origin and/or each dimension of the observations have dissimilar dispersion.\footnote{One method to at least partly ameliorate this issue would be to perform multiple iterations of the approach described by chen2016differentially. For example, preliminary bounds could be used within a DP mechanism that outputs the approximate the center of the distribution, and new DP bounds could be computed after using this first estimate to recentering the data. This process could be repeated multiple times, at the cost of requiring additional privacy-loss budget in each iteration. Note that this variant suffers from the same issue as the more basic approach; the accuracy of the final DP estimator becomes worse as the population distribution is shifted away from the origin.} Also, since the datapoints outside of $[-c,c]^d$ are removed from the dataset prior to estimation, setting $\psi$ typically requires balancing a tradeoff between bias and dispersion of the final DP estimator.
In part to bound the scope of the paper, all of the estimators proposed below use perturbation methods that result in DP estimators with log-concave distributions conditional on the data. Many alternatives that do not satisfy this property can be formulated with only minor changes to the proposed approaches, so we will point out some of these possibilities throughout the paper. However, this is also an advantageous property for a DP estimator to satisfy for two reasons. First, this condition ensures the likelihood ratio test for the population mean is monotonic, which we expect will make future work on DP inference of the proposed estimators more straightforward. Second, this condition also ensures the moments of the DP estimators exist conditional on the data. This allows for the DP estimators described here to be used with data from respondents in fairly granular geographic regions and then for summary statistics of these estimates, such as the mean or variance of these estimators across the geographic regions, to have a meaningful interpretation as unbiased estimates of their finite population counterparts. For example, in the use case described by chetty2018opportunity and using the proposed estimators, one could estimate a measure of the socioeconomic mobility in the US using a weighted mean of the estimates in each Census tract.
After outlining notation in the next subsection, the remainder of the paper is organized as follows. Section (ref) outlines the required definitions in the DP and statistical depth literature that we will use throughout the paper, and related work on DP linear regressions. Section (ref) provides tighter bounds on a parameter that is used in DP mechanisms that are based on smooth sensitivity nissim2007smooth and that use a Laplace noise distribution, which may be of independent interest to the DP community. Afterward, the proposed DP and RDP Tukey median estimators are described in Section (ref), and the proposed DP and RDP deepest regression estimators are described in Section (ref). Simulations are provided in Section (ref), and (ref) concludes.
The mechanisms described here take datasets of $n>0$ observations as input; we will denote the set of all such datasets as $\mathbb{D}^n.$ For each dataset $D\in\mathbb{D}^n,$ we will assume throughout that each respondent contributes to at most one observation in $D,$ which is assumed to be an element of $\mathbb{R}^d.$\footnote{We also assume that it makes sense to talk about the data of one respondent in isolation, which, for example, fails to hold for data on the adjacency relationships in social networks. kifer2014pufferfish provides more information on how interrelated observations impacts formal privacy guarantees.} Unless noted otherwise, we will not assume the observations in the input dataset are drawn from a population distribution. When introducing each of the formally private estimators proposed below, we will denote a candidate set, or feasible set, of estimators as $\Theta,$ and the non-private estimate as $\hat{\boldsymbol{\theta}}.$
We will use bold symbols to denote vectors, and also use $\log(x),$ where $x\in\mathbb{R}_{++},$ to denote the natural logarithm of $x.$ In addition, let $\mathbf{1}_A:\mathbb{R}^d\rightarrow \{0,1\}$ denote the indicator function of the set $A\subset\mathbb{R}^d,$ $\lfloor x \rfloor$ denote the largest integer less than or equal to $x\in\mathbb{R},$ $\lceil x \rceil$ denote the smallest integer greater than or equal to $x\in\mathbb{R},$ $x_{(k)}$ denote the $k^{\textrm{th}}$ order statistic of $\{x_i\}_{i=1}^n,$ $P_X(\cdot)$ denote the probability with respect to the random variable $X,$ $\textup{conv}(A)$ denote the convex hull of the set $A,$ and let $\textup{sign}:\mathbb{R}\rightarrow\{-1,0,1\}$ denote the sign function. Also, let $d_{\mathcal{H}}: \mathbb{D}^n \times \mathbb{D}^n \rightarrow\mathbb{Z}$ be defined so that $d_{\mathcal{H}}(D,D^\prime)$ is equal to the number of records that must be added and/or removed from $D$ to derive $D^\prime.$ Note that, since substituting one record for another using only these two operations requires adding one record and then removing one record, for any $D,D^\prime\in\mathbb{D}^n$ that differ in $k$ records, we have $d_{\mathcal{H}}(D,D^\prime)=2k.$ We will also use $\lVert \boldsymbol{x} \rVert_p$ to denote the $L^p$ norm of $\boldsymbol{x}\in\mathbb{R}^d,$ and $\lVert \boldsymbol{x} \rVert$ to denote the Euclidean norm.
The definition of an $(\epsilon,\delta)-$DP mechanism was first provided by dwork2006calibrating,dwork2006our, as described below. where all randomness is due to the mechanism
We also follow the norm in the literature and refer to $(\epsilon, 0)-$DP mechanisms as pure $\epsilon-$DP mechanisms, and we say that a given $(\epsilon, \delta)-$DP mechanism is an approximate DP mechanism when $\delta > 0.$ We will also occasionally refer to $\epsilon$ as the privacy-loss budget. Note that, since we define neighboring databases as $N(D)=\{D^\prime\in\mathbb{D}^n : d_{\mathcal{H}}(D, D^\prime) = 2 \},$ the definition above corresponds to bounded DP because the sample size of all neighbors of $D\in\mathbb{D}^n$ is fixed at $n$ and thus is bounded. The definition of unbounded DP follows from instead defining the neighboring datasets of $D\in\mathbb{D}^n$ as $N_{u}(D) = \{D^\prime\in \mathbb{D}^{n-1} \cup \mathbb{D}^{n+1} : d_{\mathcal{H}}(D, D^\prime) = 1\}.$ We use the bounded DP definition here because some of the mechanisms proposed below do not attempt to protect inferences on the sample size of the input dataset, and the norm in the DP literature is to use the bounded DP definition in these cases.
wasserman2010statistical provide an intuitive interpretation of a DP guarantee; the Neyman-Pearson lemma implies that the requirement of DP is equivalent to bounding the power of any possible hypothesis test for the null hypothesis that a given respondent's attributes are equal to a given value. Some of the many advantages of this privacy definition include invariance to post-processing, i.e., if $M:\mathbb{D}^n \rightarrow \Theta$ satisfies $(\epsilon,\delta)-$DP then so does $D \mapsto f(M(D))$ for any $f:\Theta \rightarrow \Lambda,$ and sequential composition, i.e., if $M_1:\mathbb{D}^n \rightarrow \Theta$ and $M_2:\mathbb{D}^n \rightarrow \Theta$ both satisfy $(\epsilon,\delta)-$DP then $ M(D) = (M_1(D), M_2(D))$ satisfies $(2 \epsilon, 2 \delta)-$DP. For other properties, including other forms of composition, see dwork2006calibrating, dwork2014algorithmic.
Designing an $(\epsilon,\delta)-$DP mechanism for a given function often requires a bound on its global sensitivity, which is described in the following definition, along with a related definition, which will also be used in the next subsection.
The following lemma provides a particularly simple $\epsilon-$DP mechanism known as the Laplace mechanism. Like all methods defined using global sensitivity, this mechanism requires that the global sensitivity of the function is bounded. For example, this is the case for counting queries (e.g.: a function that releases the total population within a given geographic region). However, this is not typically the case for functions that have an unbounded domain; for example, without a priori bounds on possible incomes, the following mechanism could not be used to release the average income of the residents in a geographic region. One possible way around this issue is to use a DP mechanism to approximate bounds on the observations in the sample, and then use a Laplace mechanism to release the average income of respondents with incomes that are within these bounds chen2016differentially. This approach is discussed in more detail in Section (ref).
One relaxation of $(\epsilon,\delta)-$DP is $(\epsilon, \delta,\gamma)-$random differential privacy, which is defined below. This privacy definition limits inferences between neighboring datasets that are sufficiently likely to consist of draws from the same population distribution, rather than attempting to limit inferences between all pairs of neighboring datasets. Note that this definition does not assume that the data curator knows the population distribution is in a particular class of distributions.
nissim2007smooth provide a method to decrease the sensitivity value used within a DP mechanism such as the Laplace mechanism for functions that have high global sensitivity but most often have low local sensitivity. A common example of one such function is the univariate median of values in a sample within a bounded interval, i.e., while it is possible to have a local sensitivity equal to the length of this interval, which is thus also equal to the global sensitivity, changing one record typically changes the median by a much smaller amount.
The following Lemma provides one way of deriving a $\beta-$smooth upper bound on the local sensitivity of a function, which we make use of below.
As an example, the lemma above implies that one $\beta-$smooth upper bound on the local sensitivity of the function $f(\cdot)$ can be derived by first defining the upper bound on the local sensitivity at distance $k$ as
where $\Omega(k,D)=\{f(D^\prime) : D^\prime\in\mathbb{D}^n, \; d_{\mathcal{H}}(D, D^\prime) \leq 2 k\}$ and then defining the final $\beta-$smooth upper bound as $S_f(D) = \max_{k\in \{0,1,\dots\}} \exp(-k \beta) \widetilde{A}_f^{(k)}(D).$
The following lemma provides an example of a mechanism that uses a $\beta-$smooth upper bound on local sensitivity.
tukey1975mathematics introduced the concept of halfspace depth, which can be viewed as a multivariate generalization of rank. The halfspace depth of $\boldsymbol{\theta}\in\Theta$ is the minimum number of datapoints in a halfspace with a boundary passing through $\boldsymbol{\theta},$ which is also defined below.
This measure of depth has many advantageous properties, including that it is affine equivariant, which means that, for any $A\in\mathbb{R}^{d\times d}$ with full rank and $\boldsymbol{b}\in\mathbb{R}^d,$ we have $\textup{hdepth}(\boldsymbol{\theta}, D) = \textup{hdepth}(A \boldsymbol{\theta} + b, \{A \xi + b : \xi \in D\})$ donoho1992breakdown. Also, we will make use of the following property on each halfspace-depth contour of $D,$ or the subset of $\Theta$ with half-space depth at least a given value.
Maximizing $\textup{hdepth}(\boldsymbol{\theta}, D) $ provides a measure of central tendency of a dataset that is known as the Tukey median, which is a multivariate generalization of the median. As described in the definition below, the maximizer of the halfspace depth is not generally unique, so, in cases in which this maximizer is not uniquely defined, we will define the Tukey median using a choice rule that outputs a unique element of $\Theta.$ The results in this paper only require that this choice rule outputs a point in the convex hull of data points with maximum halfspace depth, which is satisfied by any reasonable rule; one common approach in practice is to simply define the Tukey median as the arithmetic mean of the data points with maximum halfspace depth.
To describe a few advantageous properties of the Tukey median, some additional notation will be helpful. First, we will say that a dataset consists of observations in general position if there does not exist a $d-1$ dimensional hyperplane that contains more than $d$ observations. Note that, when this condition does not hold, it can be ensured with probability one by dithering the dataset, i.e., adding continuously distributed mean zero noise with a low scale to each observation; this will be discussed in more detail in Section (ref). Second, the breakdown value of an estimator is defined as the smallest possible proportion of points that must be contaminated in order to make the estimator equal to an arbitrary value, as defined below more formally.
One advantageous property of the Tukey median is that it is a robust estimator. For example, for datasets with observations in general position, its breakdown value is at least $1/(1+d),$ and one related property that we will use below is that, for any such dataset, the maximum halfspace depth is at least $\lceil n/(1+d) \rceil$ masse2002asymptotics, donoho1992breakdown. The Tukey median is also consistent at the usual parametric rate of $O(1/\sqrt{n})$ when $D$ is composed of independent observations from a population distribution. Detail on computational considerations will also be provided in the Section (ref), after outlining regression depth in the next section.
Some notation will be helpful before describing regression depth and its maximizer, the deepest regression hyperplane, which were both first proposed by rousseeuw1999regression. In this section, and all sections related to the DP deepest regression estimator below, we will suppose the dataset $D\in \mathbb{D}^n$ is composed of observations $(\boldsymbol{x_i}, y_i)\in \mathbb{R}^d.$ We also only consider regressions with an intercept throughout the paper; in other words, we suppose that each vector $\boldsymbol{\theta}\in\Theta\subset\mathbb{R}^d$ corresponds to the candidate fit $y = (1, \boldsymbol{x}^\top)\boldsymbol{\theta}.$
Before providing the general definition of regression depth, we will start by considering the case of a simple linear regression, i.e., $\Theta \subset \mathbb{R}^2.$ In this case, for any fixed value $u\notin \{x_i\}_i^n,$ the regression depth can be defined using the integers
which correspond to the number of data points either above or below the graph of the fit and either to the left or the right of $u\in \mathbb{R}^1.$ Note that rotating the graph of the fit clockwise around the point $(u, (1, u) \boldsymbol{\theta})$ to a vertical line would intersect $L^-_u(\boldsymbol{\theta}, D) + R^+_u(\boldsymbol{\theta}, D)$ observations, whereas rotating this graph counterclockwise intersects $L^+_u(\boldsymbol{\theta}, D) + R^-_u(\boldsymbol{\theta}, D)$ observations. Thus, the regression depth in this two dimensional case is given by
Regression depth can be defined in an analogous way in higher dimensional cases by replacing $u$ with a hyperplane in $\boldsymbol{x}-$space rousseeuw1999regression. We will use the definition below in this paper instead, which appears to have been first described by mizera2002depth.
This definition can be derived by considering the regression depth in the dual space. We will briefly outline this and a few related concepts, in part because the local sensitivity of the deepest regression has a straightforward geometric interpretation in the dual space. More detail on the dual space in the context of regression depth can be found in van2008efficient, rousseeuw1999regression; for interesting historical connections to related regressions, see koenker2000galton. We will use the example in Figure (ref) to introduce these concepts. The left plot in Figure (ref) depicts a set of five data points $(x_i,y_i)\inD$ in $(x,y)-$space, or the primal space. Each of these primal points can alternatively be encoded as lines in dual space using the mapping $(x_i, y_i) \mapsto \{\boldsymbol{\theta} \in \Theta : (1, x_i) \;\boldsymbol{\theta} = y_i \},$ or the set of all fits that pass through the data point in the primal space. The right plot in Figure (ref) provides the dual lines corresponding to each of the primal points. Likewise, each point $\boldsymbol{\theta} \in \Theta$ in the dual space corresponds to a single regression line in primal space using the mapping $\boldsymbol{\theta} \mapsto \{ (x, y) \in \mathbb{R}^2 : (1, x)\; \boldsymbol{\theta} = y \}.$
The act of rotating the graph of a fit to a vertical hyperplane in primal space corresponds to simply moving the point $\boldsymbol{\theta} \in \Theta$ in dual space corresponding to this fit along a ray. For example, in Figure (ref), we can see that the regression line on the left plot has a regression depth of three in two ways. First, we can count the minimum number of data points this line must intersect when rotated to a vertical line. In this case rotating this regression line counterclockwise to a vertical line around the intersection point of this line with the right hand side of the plot results in the line crossing three points in total. Second, we can alternatively count the minimum number of observation lines that a ray in dual space, originating from this regression coefficient's location, i.e., the blue point in the plot on the right side Figure (ref), must cross. The gray ray in Figure (ref) is one such example; since this ray crosses three observation lines, including the two lines at the origin of the ray, the regression depth of this fit is three.
To derive definition (ref), note that the observation line $\{ \boldsymbol{\theta} \in \Theta : y_i = (1, \boldsymbol{x_i}^\top) \boldsymbol{\theta} \}$ intersects the ray that originates from $\boldsymbol{\theta}$ and points in the direction $\boldsymbol{u} \neq \boldsymbol{0}$ if and only if there exists $r\in [0, t]$ such that
After taking the limit as $t\rightarrow \infty,$ we have
which implies Definition (ref).
Like halfspace depth, one advantageous property of regression depth is that it is affine equivariant. Also, below we will make use of the following property of regression depth contours, which are sets in dual space given by $\Omega(k, D) = \{ \boldsymbol{\theta} \in\Theta : \textup{rdepth}(\boldsymbol{\theta}, D) \geq k \}.$
One can also consider all other candidate fits to see that the fit in Figure (ref) is the fit with the highest regression depth, or the deepest regression, which is defined below. Like the Tukey median, the maximizer of regression depth need not be unique, so the definition below also uses a choice rule to ensure the deepest regression estimator is defined uniquely. Also like the Tukey median, our results only require that this choice rule outputs a final estimate in the convex hull of $\boldsymbol{\theta}\in\Theta$ with maximum regression depth. One common approach, and the approach we use in the simulations below, is to define this choice function so that it outputs the arithmetic mean of observation line intersections in dual space that have maximum regression depth.
One advantageous property of this estimator is that it is robust, with a breakdown value of at least $1/(1+d)$ mizera2002depth. In addition, the deepest regression obtains this high breakdown value without sacrificing the usual parametric rate of consistency of $O(1/\sqrt{n})$ to the conditional median of $y$ given $\boldsymbol{x};$ for more detail, see bai2008asymptotic,he1998asymptotics. Note that a similar consistency result also holds for the least absolute deviation (LAD) regression koenker1978regression, but the deepest regression has the advantage over the LAD regression of being robust to outliers in the independent variables, $\{\boldsymbol{x_i}\}_i.$ Also, one property that we will use below is that the regression depth of the deepest regression is at least $\lceil n/(1+d)\rceil,$ which, unlike the similar result described above for the halfspace depth of the Tukey median, has been shown to hold for datasets that are not in general position mizera2002depth,amenta2000regression.
In this subsection we will describe algorithms for computing halfspace and regression depth and their maximizers when $\Theta \subset \mathbb{R}^d,$ where $d\geq 2.$ Both of the notions of statistical depth can be computed at a candidate estimate $\boldsymbol{\theta}\in\Theta$ in $O(n^{d-1} \log(n))$ time, and faster approximations for these methods are also available langerman2000optimal,rousseeuw1998computing.
In the case of the Tukey median, simply computing the halfspace depth at each datapoint provides an algorithm for computing the Tukey median in $O(n^d\log(n))$ time. Faster approaches are also available when $d=2$ miller2001fast,langerman2003optimization. struyf2000high also provide an approximation method for the Tukey median when the sample size and/or the dimension is prohibitively large.
In the case of regression depth, there are $O(n^d)$ candidate fits that pass through $d$ data points to consider, so the na\"ive approach of computing the regression depth for each of these fits has a time complexity of $O(n^{2 d - 1} \log(n)),$ and is most often prohibitively computationally expensive when $d>2.$ langerman2003complexity provide a method for approximating the deepest regression that runs in $O(n\log(n))$ time when $d=2.$ van2002deepest provide a method for approximating to the deepest regression in use cases in which the sample size and/or the dimension is higher.
We emphasize that the first four differentially private mechanisms provided below actually require more than the maximizer of a notion of statistical depth; in order to also compute a $\beta-$smooth upper bound on the local sensitivity of the estimator, these mechanisms also require upper bounds on diameters of the depth contours defined in Lemma (ref) and Lemma (ref).
Computationally efficient methods to compute upper bounds on these values are left as a topic for future research. For the purposes of these first four differentially private methods in this paper we simply use the two na\"ive approaches in order to explicitly find these diameters. Since this implies that the deepest regression mechanisms provided below are primarily applicable when $d=2,$ Section (ref) provides a mechanism for the more computationally efficient approximation of the deepest regression provided by van2002deepest.
This section mainly focuses on providing some pointers into the DP literature for linear regression and Tukey median estimators. The only paper that we are of that provides a DP Tukey median estimator is ramsay2021differentially, which uses a formulation based on an exponential mechanism rather than on the $\beta-$smooth sensitivity of the Tukey median, as is done below. One such class of DP estimators are those that pass sufficient statistics of a linear regression through a DP primitive mechanism; see for example, foulds2016theory,mcsherry2009differentially,vu2009differential,wang2018revisiting,dwork2014analyze. A related class of DP estimators include the methods proposed by dwork2009differential,alabi2020differentially, which are based on the Theil-Sen simple linear regression estimator theil1950rank,sen1968estimates. bassily2014private also propose a DP OLS method that estimates the regression parameters using a DP gradient descent method. One line of work that can be used to formulate a particularly broad class of DP estimators is based on perturbing the objective function of a convex optimization problem kifer2012private,awan2021structure. There has also been work on DP inference; see for example, barrientos2019differentially,sheffet2017differentially,awan2021structure,karwa2018finite.
Since some of our results are related to the input requirement for most DP estimators of a priori bounds on the observations, it is worth highlighting a few papers related to this issue. First, as described above, chen2016differentially describe a DP mechanism that can be used to estimate a bounding box of the form $[-c,c]^d \subset \mathbb{R}^d$ that contains a fixed proportion of the datapoints. These bounds can be used as input to subsequent DP mechanisms after removing all observations that are outside of this bounding box from the dataset. More detail on this approach is also given in the next subsection and Section (ref). Second, karwa2018finite provide DP estimators for the univariate mean of the data without the requirement of input bounds on the observations. While the $(\epsilon,\delta,\gamma)-$RDP estimators proposed below also do not require input bounds on the data or the estimator, the $(\epsilon,\delta,\gamma)-$RDP estimators proposed below are different from the (pure and approximate) DP estimators proposed by karwa2018finite for a few reasons. First, the accuracy of the approach described by karwa2018finite is dependent on whether or not the input dataset satisfies their assumption that it is composed of independent draws from a Gaussian distribution; in contrast, the $(\epsilon,\delta,\gamma)-$RDP estimators proposed below do not require distributional assumptions on the population distribution. Second, since RDP is a relaxation of the privacy guarantee provided by DP, the privacy guarantee provided by the mechanisms proposed by karwa2018finite is stronger than the corresponding guarantee provided by the estimators proposed below. Third, in the case of the RDP estimators proposed below, our avoidance of distributional assumptions requires an assumption on the sample size being above a given threshold, which is not required by the approach described by karwa2018finite. Later work also generalized this approach to the multidimensional case, while maintaining the normality assumption on the population distribution; see for example, kamath2022private.
This section provides improved bounds on the maximum value of $\beta$ that can be used in a smooth-sensitivity-based mechanism with noise drawn from a Laplace distribution while still satisfying $(\epsilon, \delta)-$DP. The corresponding bounds on $\beta$ that are provided by nissim2007smooth in this case, which are also provided in Lemma (ref), prioritize tractability over tightness. For example, when $d=2$ and $\delta \leq 1/10^6,$ the bound proposed by nissim2007smooth for this case is slightly less than half of the value provided here. The result that provides the improved lower bound on $\beta,$ i.e., Theorem (ref), follows from the next theorem, which provides three upper bounds on a certain quantile of the sum of exponential random variables, given in ((ref))-((ref)). In the simulations provided in Section (ref), bound ((ref)) is used to set $\beta$ to the tightest bound we derive here, as described in Theorem (ref).
The next result describes how Theorem (ref) can be used to derive possible values for $\beta.$
As described previously, in the simulations presented below, we use the tightest bound implied by this result, which is $\beta = \epsilon/(2 q(1-\delta; d, 1)).$ The next two sections will provide the DP mechanisms that make use of the $\beta-$smooth upper bounds on local sensitivity.
The next two subsections will introduce the proposed formally private location and regression estimators, respectively. Specifically, both sections will start by introducing an $(\epsilon,\delta)-$DP estimator that, as described above, requires the user to input an a priori feasible set for the estimator, and the proposed $(\epsilon,\delta,\gamma)-$RDP mechanisms are based on estimating a suitable feasible set.
Our results on the $(\epsilon,\delta,\gamma)-$RDP estimators will also provide sufficient conditions for these mechanisms to satisfy information-set-invariant nonrevealing (ISIN). This property ensures that, even for cases in which a privacy attacker knows all but one observation, the mechanism does not reveal any unknown observations with certainty, as defined in more detail below.
Note that all pure $\epsilon-$DP mechanisms satisfy this property. In practice, most commonly used approximate $(\epsilon,\delta)-$DP mechanisms satisfy this property as well, but it is also possible to construct counterexamples that do not. For example, consider the $(\epsilon,\delta)-$DP mechanism that outputs approximate total populations of each Census block. Suppose the true counts are given by $\boldsymbol{\theta}\in\mathbb{N}^d_+.$ Note that the global sensitivity of this query is $2,$ so one $(\epsilon,\delta)-$DP mechanism could be defined using $\boldsymbol{z}\in\mathbb{R}^d,$ where $\boldsymbol{z}[i]\sim \textrm{Laplace}(0, 2/\epsilon),$ and $b>0$ defined so that $P(\lVert \boldsymbol{z} \rVert_1 \leq b) = 1-\delta.$ The output of the mechanism can be defined as $\boldsymbol{\theta} + \boldsymbol{z}$ when $\lVert \boldsymbol{z} \rVert_1 \leq b$ and $\boldsymbol{\theta}$ directly otherwise. Since the Laplace distribution is absolutely continuous, whenever this mechanism's output is a vector of integers, a privacy attacker would know that this output is the vector of true counts with probability one.
Despite the fact that there exists examples of $(\epsilon,\delta)-$DP mechanisms that do not satisfy ISIN, these mechanisms are often avoided in practice. However, as described above, the original example of an $(\epsilon,\delta,\gamma)-$RDP estimator provided by hall2013random, which was a sparse histogram estimator, did not satisfy ISIN. This is because the input dataset itself was used to estimate the sparsity pattern of the population distribution, and no noise was added to histogram counts that were very likely to be zero. In cases in which a privacy attacker knows the data of all but one respondent and the unknown respondent changes the sparsity pattern of the true histogram, the data of the unknown respondent is revealed with certainty. In contrast, ensuring the $(\epsilon,\delta,\gamma)-$RDP estimators proposed here satisfy this property is straightforward because we use the input dataset for a more limited purpose.
In this section, we describe both an $(\epsilon,\delta)-$DP and an $(\epsilon,\delta,\gamma)-$RDP Tukey median estimators. Afterward, to take advantage of tighter probability bounds in the univariate case, we also provide a separate $(\epsilon,\delta,\gamma)-$RDP median estimator. Note that we provide more detail in the proofs in this section than the next section, which provides the corresponding results for the $(\epsilon,\delta)-$DP and $(\epsilon,\delta,\gamma)-$RDP deepest regression estimators, since the results for both of these estimators follow from identical arguments.
Even though the results on the deepest regression mechanisms follow from the same arguments as the corresponding results for the Tukey median mechanisms, the assumptions are not identical. Specifically, in the case of our proposed $(\epsilon,\delta)-$DP and $(\epsilon,\delta,\gamma)-$RDP Tukey median estimators when $d\geq 2,$ we will require the additional assumption that the dataset consists of datapoints in general position. This can be ensured without assuming the input dataset consists of draws from an absolutely continuous population distribution using dithering, as is also described in Section (ref). Specifically, each $\boldsymbol{x_i}\inD$ can be redefined as $\tilde{\boldsymbol{x_i}} = \boldsymbol{x_i} + \boldsymbol{z_i},$ where $\boldsymbol{z_i}\in\mathbb{R}^d$ is a mean zero random variable with a distribution that has a low dispersion and is absolutely continuous. Note that the result of using this technique is that we are providing privacy guarantees to the records given by $\{\tilde{\boldsymbol{x_i}}\}_i,$ for a given set of realizations of the noise random variable $\{\boldsymbol{z_i}\}_i,$ rather than directly to $\{\boldsymbol{x_i}\}_i.$ The benefit of this approach is that it allows us to constrain our view of the universe of input datasets to those that are in general position; however, in use cases that require using sequential composition to release the output of several formally private mechanisms, including a formally private Tukey median estimate, ensuring the privacy guarantee provided by each mechanism holds on the same input dataset would require using $\{\tilde{\boldsymbol{x_i}}\}_i$ as the input dataset in all of these mechanisms. For this reason, after the next Lemma we describe how this assumption can be removed entirely using an alternative approach. In some cases, particularly when the a priori bounds on $\Theta$ are relatively tight and/or the sample size is sufficiently large, the negative impact of using this alternative on the accuracy of the final DP estimator is likely small or even entirely nonexistant.
Both of the mechanisms provided in this section will use the following $\beta-$smooth upper bound on the local sensitivity of the Tukey median estimator.
Algorithm (ref) provides a mechanism to approximate the Tukey median that satisfies $(\epsilon,\delta)-$DP, which is an immediate consequence of the previous lemma and the main result of nissim2007smooth, as stated in the following corollary. Note that it is also possible to use the $\beta-$smooth upper bound on the local sensitivity of the Tukey median provided in Lemma (ref) to formulate pure $\epsilon-$DP mechanisms, at the cost of the distribution of the resulting estimator not having exponential tails, as is described by nissim2007smooth in more detail.
The next theorem also provides a method to define the set $\Theta$ that is the input into DPTukeyMedian$(\cdot)$ in a data dependent manner while satisfying $(\epsilon, \delta, \gamma)-$RDP when $d\geq 2.$ Note that we treat the case in which $d=1$ separately in Theorem (ref).
The final result of this section improves on Theorem (ref) for the case in which $d=1.$ Note that the assumption that the datapoints are in general position is not required for this Theorem.
In this section, we will provide deepest regression mechanisms that satisfy formal privacy guarantees. We will start by providing the regression depth counterpart to Lemma (ref).
Algorithm (ref) provides an $(\epsilon,\delta)-$DP mechanism for approximating the deepest regression, which is an immediate consequence of the Lemma (ref) and the main result of nissim2007smooth, as stated in the following corollary. Note that the same comment preceding Corollary (ref) also holds in this case; it is also straightforward to use Lemma (ref) to formulate pure $\epsilon-$DP deepest regression estimators nissim2007smooth.
The next theorem also provides a method to set the input $ \Theta$ of DPDeepestReg$(D, \Theta, \epsilon, \delta)$ using the data directly while satisfying $(\epsilon, \delta, \gamma)-$RDP. This mechanism can be made to satisfy ISIN by simply expanding this feasible set input, as described in more detail in the theorem. The simulations in Section (ref) provide evidence that even a large expansion of this feasible set input may not have any impact on the distribution of the resulting estimator when the sample size is sufficiently large.
As described previously, the deepest regression methods described in the previous section are most often limited to the case in which $d=2,$ so this section provides a computationally efficient method for computing an $(\epsilon,\delta)-$DP approximate deepest regression estimate called Medsweep that was proposed by van2002deepest. In this section it will be helpful to define $X\in\mathbb{R}^{n\times (d-1)}$ so that $X[i,\cdot]=\boldsymbol{x_i}$ and $\boldsymbol{y}\in \mathbb{R}^n$ so that $\boldsymbol{y}[i]=y_i.$
Informally, the method can be viewed as similar to applying the Frish-Waugh-Lovell theorem when finding an ordinary least squares estimator, in the sense that the method proceeds by parsing (or “sweeping") $X[\cdot,i]$ out of $X[\cdot, k]$ by approximating the univariate regression of $X[\cdot,i]$ on $X[\cdot, k]$ and then redefining $X[\cdot,k]$ as the residual of this regression. Afterward, a similar approach is used to iteratively parse each column $X[\cdot,k]$ out of $\boldsymbol{y}.$
There are a few possible ways to formulate DP Medsweep mechanisms. The simplest approach is to replace each median evaluation within Medsweep with a DP median mechanism. Here we will describe a slightly more involved formulation that has the advantage of providing a DP estimate with exponential tails conditional on the data. Specifically, this requires formulating DP variants of two methods that are used within Medsweep. The first is simply the median; below we use the $\beta-$smooth sensitivity of the median provided by nissim2007smooth to define the scale of a Laplace mechanism, as described in equation ((ref)).
The second method is $r:\mathbb{R}^n\times \mathbb{R}^n \rightarrow \mathbb{R},$ and is defined as
The following lemma provides the $\beta-$smooth upper bound on the local sensitivity of $r(\cdot)$ when its output is truncated to the range $[L, U]\subset\mathbb{R}.$ This method is used in two places in Medsweep algorithm, with inputs $\boldsymbol{u},\boldsymbol{v}\in \mathbb{R}^n$ are given by either columns of $X$ or one column of $X$ and the vector $\boldsymbol{y}.$ Note that the final DP mechanism actually imposes bounds on the output of these intermediate mechanisms, rather than the final estimator.
The function $\widetilde{A}^{(k)}_r(\boldsymbol{u}, \boldsymbol{v})$ can alternatively be defined using a concept from computational geometry literature; some additional notation will be helpful to describe this connection. We will let $\Gamma$ denote a set of curves in $\mathbb{R}^3,$ with curve $i\in\{1,\dots,n\}$ defined as $ \gamma(i, \tilde{u}, \tilde{v}) = (\boldsymbol{u}[i]-\tilde{u} )/(\boldsymbol{v}[i]-\tilde{v}).$ Let $\mathcal{A}(\Gamma)$ denote the arrangement of $\Gamma,$ which is defined as the collection of graphs of functions in $\Gamma.$ More detail on arrangements and the other concepts from the computational geometry that we briefly describe here can be found in edelsbrunner1987algorithms. We will define the level of a point $p\in\mathbb{R}^3$ as the number of curves in $\Gamma$ that either pass below or through $p,$ and define the $k-$level of $\mathcal{A}(\Gamma),$ or $L(k, \Gamma),$ as the points on one or more curves in $\Gamma$ with level $k.$ Also, a $\zeta-$approximate $k-$level of an arrangement of curves is defined as a curve that is within the $k-\zeta$ and $k+\zeta$ levels of an arrangement, which we will denote by $L_\zeta(k, \Gamma).$ Using this notation, we have,
Unfortunately, we are not aware of a computationally efficient method of computing approximations of all levels of $\{L(k, \Gamma)\}_k$ for this use case, in which the curves in $\Gamma$ are nonlinear and defined in $\mathbb{R}^3.$\footnote{agarwal1990partitioning provides a related method, which provides approximations to all levels when $\Gamma$ is a collection of lines in $\mathbb{R}^2.$ This can used within a method to compute a $\beta-$smooth upper bound on the local sensitivity using some of the same techniques in the approach described below. For example, upper bounds on the order statistics of $\{\gamma(i, \tilde{u}, \tilde{v})\}_{i=1}^n$ can be computed for $\tilde{v}$ in each partition element $q_j = (j c, (j+1) c],$ where $c>0,$ by constructing lines that are upper bounds of $\tilde{u} \mapsto \gamma(i,\tilde{u},\tilde{v})$ over all $\tilde{v}\in q_j.$ However, there are some issues with using this approach in practice, such as the large constants in the asymptotic time complexity of the approach described by agarwal1990partitioning.} Algorithm (ref) provides an alternative approach, which is also used in the simulations provided in the next section. This approach uses the choice parameters $c_u,c_v>0$ to define the partitions of $\mathbb{R}^1$ given by $q_u=\{\cdots, [c_u (j-1), c_u j), [c_u j, c_u (j + 1)), \cdots\}_{j=-\infty}^\infty$ and $ q_v=\{\cdots, [c_v (j-1), c_v j), [c_v j, c_v (j + 1)), \cdots\}_{j=-\infty}^\infty.$ After defining a partition of $[\boldsymbol{u}_{(m-k)}, \boldsymbol{u}_{(m+k)}] \times [\boldsymbol{v}_{(m-k)}, \boldsymbol{v}_{(m+k)}] \subset \mathbb{R}^2$ as
we bound each of the order statistics of $\{\gamma(i,\tilde{u}, \tilde{v})\}_i,$ for each $(\tilde{u}, \tilde{v})$ in the partition element $ c \in B(k,\boldsymbol{u},\boldsymbol{v}),$ between the corresponding order statistic of $ \{ \inf_{(\tilde{u}, \tilde{v}) \in c} \gamma(i,\tilde{u}, \tilde{v}) \}_i$ and $\{\sup_{(\tilde{u}, \tilde{v}) \in c} \gamma(i,\tilde{u}, \tilde{v})\}_i.$ However, if $c_u,c_v$ are defined to be $O(1/n^{\alpha })$ for a choice parameter $\alpha > 0$ and $\{(\boldsymbol{u}[i], \boldsymbol{v}[i])\}_i$ consists of independent and identically distributed observations, this algorithm's output still converges to zero and has a time complexity of $O_p(n^{1 + 2\alpha } \log(n)),$ as described in the next Lemma. In the simulations in the next section we defined $c_u$ and $c_v$ as $8/n^{3/4},$ which results in a time complexity of $O_p(n^{2.5}\log(n)).$
Algorithm (ref) provides pseudocode for the approximate DP variant of Medsweep; the following theorem provides the privacy guarantee of this mechanism.
The primary advantage of the approach used in Algorithm (ref) is that the noise added to the output has exponential tails. However, there are several alternatives to this approach that are worth mentioning as well, particularly in cases in which noise with exponential tails is not required. First, as described above, one can alternatively formulate a DP variant of Medsweep by replacing each median evaluation with a DP median mechanism. This can be done using a variety of DP median mechanisms, including those that use smooth sensitivity, as described by nissim2007smooth, or an approach based on the exponential mechanism provided by mcsherry2007mechanism. Second, since the mechanism below uses several calls to primitive mechanisms, especially in the case of multivariate regressions, $\rho-$concentrated differential privacy ($\rho-$CDP) variants are also worth mentioning because this privacy guarantee provides tight sequential composition. This privacy guarantee was first described by bun2016concentrated; also, bun2019average describe how $\rho-$CDP mechanisms can be formulated using $\beta-$smooth upper bounds on local sensitivity.
This section provides simulations to compare the accuracy of DPDeepestReg$(\cdot)$ and DPMedsweep$(\cdot)$ with three other DP regression methods for three different data generating processes (DGPs). Specifically, for each $F\in \{N(\mu=0, \sigma^2=1), \textup{Laplace}(\mu=0, b=1), t(\mu=0, \nu=3)\},$ each $i\in\{1,\dots,n\},$ and each $n\in \{50,100,200,300,400,500\},$ we define the dataset $D$ as $\{(x_i, y_i)\}_{i=1}^n,$ where
Note that, since each $F\in \{N(\mu=0, \sigma^2=1), \textup{Laplace}(\mu=0, b=1), t(\mu=0, \nu=3)\}$ is symmetric, the mean of $y$ conditional on $x$ is equal to the median of $y$ conditional on $x,$ so the DP conditional mean and median estimators described in the next subsection estimate the same population value. These simulations were performed for each $\epsilon\in \{4,12\},$ and, for each $(\epsilon,\delta)-$DP mechanism, we set $\delta$ to $10^{-6}.$
For each simulation iteration, we estimate $\boldsymbol{\theta}$ using DPDeepestReg$(\cdot),$ DPMedsweep$(\cdot),$ and the three classic DP estimators described in the next subsection. As described previously, our aim with the implementation choices described in the next section was to err on the side of understating the relative accuracy of the two proposed estimators, particularly for choices related to both datapoint bounds (in the case of two of the three classic DP methods) and estimator bounds (in the case of the proposed approaches), which are described in more detail in the final paragraph of the next subsection.
The implementation for DPMedsweep$(\cdot)$ used the pseudocode given in Algorithm (ref). The implementation of DPDeepestReg$(\cdot)$ used the pseudocode given in Algorithm (ref); to compute the deepest regression estimator, we simply computed the depth of all candidate regressions that pass through each pair of data points, using the function rdepth2$(\cdot)$ from segaert2017mrfdepth that was rewritten in python.
All three of the existing DP regression estimators are based on python code from alabi2020differentially, which is available at \url{https://github.com/anonymous-conf/dplr}. The first method is DPMedTS_exp$(\cdot),$ which is a differentially private Theil-Sen estimator theil1950rank,sen1968estimates. This $\epsilon-$DP mechanism is based on the mechanism provided by dwork2009differential, but also incorporates subsequent results in the DP literature that improve accuracy of $\epsilon-$DP median mechanisms.
Second, NoisyStats$(\cdot)$ is an $\epsilon-$DP OLS estimator. This mechanism works by passing sufficient statistics for the OLS estimator through a Laplace mechanism. One disadvantage of this approach is that, whenever the $\epsilon-$DP approximation of var$(x) = \sum_i (x_i- \overline{x})^2/n$ is negative, this mechanism uses $2/3$ of the total privacy-loss budget without actually providing an estimate. For the purposes of calculating the error metrics that are reported below, these simulation iterations were simply removed before computing these error metrics; the count of the simulation iterations for which no output was returned from this mechanism is provided in Table (ref). Note that this occurred most often for smaller sample sizes, the $\epsilon=4$ simulations, and the DGPs with thicker tails that we consider.
Third, DPGradDesc$(\cdot)$ is an $(\epsilon,\delta)-$DP mechanism that provides a DP OLS estimate using a DP gradient descent method. Like the methods proposed here, we used $\delta=1/10^6$ for this mechanism. We initialized these gradient descent steps at the origin and chose an iteration limit of 30. and we also set the iteration limit to 30. Another choice parameter was set so that the gradient was clipped to the unit cube $[-10, 10]^2.$ These choices appeared to increase accuracy in our simulations relative to their default values, but it is possible that other choices would provide a further improvement.
These three existing mechanisms were designed to return estimates of the dependent variable at two values of the independent variable, so we also modified DPDeepestReg$(\cdot)$ and DPMedsweep$(\cdot)$ to output estimates of the dependent variable at two values of the independent variable. Specifically, in each simulation iteration, we saved the errors of the dependent variable estimate at the 0.25 and 0.75 quantiles of the random variable $x_i$ for each of the five DP estimators considered, and, after 100 simulation iterations, we computed the mean squared error and median absolute error for each of these estimators, which are provided in Figures (ref) and (ref) respectively. These results are also discussed in more detail in the next subsection.
The first two existing DP mechanisms described above require bounds on $\{(x_i,y_i)\}_i.$ To define these bounds, we use an approach that was inspired by the DP mechanism provided by chen2016differentially, which, as described above, outputs an approximation of the bounding box given by $[-c, c]^2,$ where $c\in \mathbb{R}_+$ is as small as possible such that at least a proportion $\psi\in (0,1)$ of the data points are in $[-c, c]^2.$ Setting this choice parameter in practice requires balancing two tradeoffs. First, higher values typically result in larger datapoint bounds, which leads to a higher global sensitivity and thus a higher variance of the noise in these mechanisms. Second, since all datapoints outside of these bounds are ignored in these mechanisms, lower values tend to decrease the statistical efficiency and increase the bias of the DP estimator. Rather than use this DP mechanism, we simply set the bounds on $\{(x_i,y_i)\}_i$ using this definition of $[-c, c]^2$ directly, with $\psi=0.98$ in the simulation results presented in this section. Since this choice parameter did have a significant impact on the performance of both of these classical DP estimators, Appendix (ref) provides the simulation results for both $\psi=0.95$ and $ \psi=1.$ The approach taken here is intended to simulate an ideal realization of the output of the preliminary DP mechanism provided by chen2016differentially. Note that we do not account for the privacy-loss budget that would be required to compute these data-dependent bounds in an actual use case in order to err on the side of overstating the accuracy of these existing methods for a given total privacy-loss budget. In contrast, for the two proposed DP estimators, we use the non-data-dependent bounds on $\boldsymbol{\theta}$ of $\Theta=[-50, 50]^2.$ The impact of this choice is described in the final paragraph of the next subsection. As described in the paragraph preceding Lemma (ref) in more detail, we defined $c_u,c_v$ as $8/n^{0.75},$ which results in DPMedsweep$(\cdot)$ having a time complexity of $O(n^{2.5}\log(n)).$
Before discussing the simulation results in detail, it may be worth pointing out our reasons for including both the mean squared error and the median absolute error of these methods. Both of these error metrics are examples of loss functions that one may seek to minimize with the choice of an estimation strategies. If this is done, the mean squared error may be regarded as the more prudent choice, in the sense that the upper tail of the absolute errors are given more weight. However, since NoisyStats$(\cdot)$ has an infinite expected absolute error, the mean squared error estimates for this mechanism do not converge in the number of simulation iterations. Thus, we choose to include median absolute error metrics solely because the population counterparts of these metrics exist for each mechanism, rather than because the invariance of median absolute error on the magnitude of the largest errors is an advantage in the context of selecting an estimator.
In the case of both mean squared error and median absolute error, one interesting feature of the simulation results is the consistency of DPGradDescent$(\cdot)$ across DGPs, values of $\epsilon,$ and sample sizes. This mechanism does not appear to exhibit a gain in accuracy as the sample sized increases; although, as described previously, this may be due to the fact that we only performed limited tuning of the choice parameters of this method. This mechanism also appears to outperform other estimators in terms of median absolute error and mean squared error in all DGPs when $\epsilon=4$ and $n\in\{50,100\}.$ In contrast, in the case of DPMedTS_exp$(\cdot),$ both the error metrics we consider are more obviously decreasing in the sample sizes for most of the DGPs and $\epsilon$ combinations we consider. It is interesting to note that even the exponential tails of the Laplace distribution appear to be sufficiently thick to make the impact of sample size less apparent than for the normal distribution for DPMedTS_exp$(\cdot).$ Since the NoisyStats$(\cdot)$ estimator does not have finite expected squared error, we will instead focus on the median absolute errors for this estimator. The impact of the tail behavior of the DGP on the statistical efficiency, i.e., the rate at which errors decrease as the sample size increases, appears to be less obvious in this case. For both DPMedTS_exp$(\cdot)$ and NoisyStats$(\cdot),$ higher values of $\epsilon$ also appear to increase the statistical efficiency.
The two estimators proposed in this paper appear to provide similar performance to one another, but DPDeepestReg$(\cdot)$ appears to provide a slight gain in accuracy over DPMedsweep$(\cdot).$ These two estimators also appear to compare favorably to the classic DP regression estimators in these simulations when either the sample size or epsilon is sufficiently large. Also, note that this accuracy improvement is most notable in cases in which the tails of the DGP are thicker.
To assess computational costs of the proposed approaches we also measured the median runtime for each sample size. On an ordinary laptop, in each simulation iteration with a sample size of 500, the median runtime was 5.8 seconds for DPDeepestReg$(\cdot)$ and 1.9 seconds for DPMedsweep$(\cdot).$ Since the vast majority of the runtime of DPMedsweep$(\cdot)$ was spent computing the smooth sensitivity, the runtime of this mechanism can be improved by reducing $\alpha,$ as defined in in Lemma (ref), below $3/4,$ at the cost of increasing the noise scale used within this mechanism.
One advantageous feature of the proposed approaches in these simulations is that $\Theta$ can be defined as progressively larger sets as the sample size diverges without changing the distribution of the DP estimators. In other words, when the sample size is sufficiently large, the bound on the local sensitivity at distance $k$ multiplied by $\exp(-k \beta),$ i.e., $\exp(-k \beta) \widetilde{A}^{(k)}_{\hat{\boldsymbol{\theta}}}(D),$ typically obtains its maximum value at a relatively low value of $k,$ while the choice of $\Theta$ only impacts the value of $\exp(-k \beta) \widetilde{A}^{(k)}_{\hat{\boldsymbol{\theta}}}(D)$ when $k$ is large. In contrast, for the three classical approaches, the choice of bounds has a significant impact on the distribution of the output in all sample sizes. To explore this feature in the context of these simulations, we computed the minimum required diameter of $\Theta,$ defined as $\textup{diam}(\Theta) = \max_{\boldsymbol{\theta}_1,\boldsymbol{\theta}_2\in\Theta} \lVert \boldsymbol{\theta}_1-\boldsymbol{\theta}_2 \rVert_1,$ that would be required to change the scale of the Laplace noise used within DPDeepestReg$(\cdot)$ for at least one of the simulation iterations of any of the three DGPs we consider in these simulations, for each sample size and each $\epsilon.$ As shown in Table (ref), the resulting minimum diameter values appear to be increasing at an exponential rate in $n$ in these simulations, with a higher rate when $\epsilon=12.$
This paper proposes several approximate DP and RDP estimators based on maximizing notions of statistical depth, including the Tukey median and the deepest regression estimator. An approximate DP mechanism for the Medsweep regression estimator is also proposed. Afterward, simulations are provided for both of the proposed approximate DP regression estimators, and these estimators appear to compare favorably to classic DP regression estimators when either the sample size or epsilon is large.
One reason the input requirement of bounds on the data or the estimator itself is an area of emphasis in this paper is that we are not aware of approximate DP location or regression estimators that provide a reasonable balance between accuracy and privacy in the absence of strong assumptions on the population distribution, as described in Sections (ref) and (ref). These assumptions are unlikely to hold in practice for many use cases of interest, especially for use cases in which one or more variables exhibit thick tail behavior, and we are also not aware of a way to check these assumptions using a DP mechanism without making other strong distributional assumptions.
There are also a number of areas for further research in the intersection of DP estimators and statistical depth. First, rousseeuw1999regression provide a quantile regression estimator based on regression depth, along the lines of koenker1978regression, and it appears that similar techniques to the ones used here can be used to formulate $\beta-$smooth upper bounds on the local sensitivity of these estimators. Second, improving the bound described in Theorem (ref) would also be advantageous for multivariate regression estimation. Third, since inference is often one of the main goals when using regressions, one important area for further work is the development of methods related to DP inference on the estimators described here.
The author thanks John M. Abowd, Brian Finley, Daniel Kifer, Roger Koenker, and Tucker McElroy for their insightful comments.