How to Find a Point in the Convex Hull Privately

Several papers studied the PIP problem for $d=1$. In particular, three different constructions with sample complexity $2^{O(\log^{*}|G|)}$ were presented in [3, 4, 5] (for $d=1$). Recently, Kaplan et al. [11] presented a new construction with sample complexity $O((\log^{*}|G|)^{1.5})$ (again, for $d=1$). Bun et al. [5] gave a lower bound showing that every differentially private algorithm for this task must have sample complexity $\Omega(\log^{*}|G|)$. Beimel et al. [2] incorporated a dependency in $d$ to this lower bound, and showed that every differentially private algorithm for the PIP problem must use at least $n=\Omega(d+\log^{*}|G|)$ input points.

For the case of pure differential privacy (i.e., $\delta=0$), a lower bound of $n=\Omega(\log X)$ on the sample complexity follows from the results of Beimel et al. [1]. This lower bound is tight, as an algorithm with sample complexity $n=O(\log X)$ (for $d=1$) can be obtained using a generic tool in the literature of differential privacy, called the exponential mechanism [15]. We sketch this application of the exponential mechanism here. (For a precise presentation of the exponential mechanism see Section 2.1.) Let $G=\{0,\frac{1}{X},\frac{2}{X},\dots,1\}$ be our (1-dimensional) grid within the interval $[0,1]$, and let $S$ be a multiset containing $n$ points from $G$. The algorithm is as follows.

1. 1.

   For every $y\in G$ define the score $q_{S}(y)=\min\left\{|\{x\in S\mid x\geq y\}|,\;|\{x\in S\mid x\leq y\}|\right\}.$

2. 2.

   Output $y\in G$ with probability proportional to $e^{{\varepsilon}\cdot q_{S}(y)}$.

Intuitively, this algorithm satisfies differential privacy because changing one element of $S$ changes the score $q_{S}(y)$ by at most $\pm 1$, and thus changes the probabilities with which we sample elements by roughly an $e^{{\varepsilon}}$ factor. As for the utility analysis, observe that $\exists y\in G$ with $q_{S}(y)\geq\frac{n}{2}$, and the probability of picking this point is (at least) proportional to $e^{{\varepsilon}n/2}$. As this probability increases exponentially with $n$, by setting $n$ to be big enough we can ensure that points $y^{\prime}$ outside of the convex hull (those with $q_{S}(y^{\prime})=0$) get picked with very low probability.

Beimel et al. [2] observed that this algorithm extends to higher dimensions by replacing $q_{S}(y)$ with the Tukey depth ${\sf td}_{S}(y)$ of the point $y$ with respect to the input set $S$ (the Tukey depth of a point $y$ is the minimal number of points that need to be removed from $S$ to ensure that $y$ is not in the convex hull of the remaining input points). However, even though there must exist a point $y\in{\mathbb{R}}^{d}$ with high Tukey depth (at least $n/(d+1)$; see [14]), the finite grid $G\subseteq{\mathbb{R}}^{d}$ might fail to contain such a point. Hence, Beimel et al. [2] first refined the grid $G$ into a grid $G^{\prime}$ that contains a point with high Tukey depth, and then randomly picked a point $y$ from $G^{\prime}$ with probability proportional to $e^{{\varepsilon}\cdot{\sf td}_{S}(y)}$. To compute the probabilities with which grid points are sampled, the algorithm in [2] explicitly computes the Tukey depth of every point in $G^{\prime}$, which, because of the way in which $G^{\prime}$ is defined, results in running time of at least $\Omega(|G|^{d^{2}})=\Omega(X^{d^{3}})$ and sample complexity $n=O(d^{3}\log|G|)=O(d^{4}\log X)$. Beimel et al. then presented an improvement of this algorithm with reduced sample complexity of $n=O(d^{2.5}\cdot 8^{\log^{*}|G|})$, but the running time remained $\Omega(|G|^{d^{2}})=\Omega(X^{d^{3}})$.

We give an approximate differentially private algorithm for the private interior point problem that runs in $O(n^{d})$ time,¹¹1When we use $O$-notation for time complexity we hide logarithmic factors in $X$, $1/{\varepsilon}$, $1/\delta$, and polynomial factors in $d$. We assume operations on real numbers in $O(1)$ time (the so-called real RAM model). and succeeds with high probability when the size of its input is $\Omega(\frac{d^{4}}{{\varepsilon}}\log\frac{X}{\delta})$. Our algorithm is obtained by carefully implementing the exponential mechanism and reducing its running time from $\Omega(|G|^{d^{2}})=\Omega(X^{d^{3}})$ to $O(n^{d})$. We now give an informal overview of this result.

To avoid the need to extend the grid and to calculate the Tukey depth of each point in the extended grid, we sample our output directly from $[0,1]^{d}$. To compute the probabilities with which we sample a point from $[0,1]^{d}$ we compute, for each $k$ in an appropriate range, the volume of the Tukey region of depth $k$, which we denote as $D_{k}$. (That is, $D_{k}$ is the region in $[0,1]^{d}$ containing all points with Tukey depth exactly $k$.) We then sample a value $k\in[n]$ with probability proportional to ${\rm Vol}(D_{k})\cdot e^{{\varepsilon}k}$, and then sample a random point uniformly from $D_{k}$.

Observe that this strategy picks a point with Tukey depth $k$ with probability proportional to ${\rm Vol}(D_{k})\cdot e^{{\varepsilon}k}$. Hence, if for a “large enough” value of $k$ (say $k\geq\frac{n}{cd}$ for a suitable absolute constant $c>1$) we have that ${\rm Vol}(D_{k})$ is “large enough”, then a point with Tukey depth $k$ is picked with high probability. However, if ${\rm Vol}(D_{k})=0$ (or too small) then a point with Tukey depth $k$ is picked with probability zero (or with too small a probability). Therefore, to apply this strategy, we derive a lower bound on the volume of every non-degenerate Tukey region, showing that if the volume is non-zero, then it is at least $\Omega\left(1/X^{d^{3}}\right)$.

There are two issues here. The first issue is that the best bound we know on the complexity of a Tukey region is $O(n^{(d-1)\lfloor d/2\rfloor})$, so we cannot compute these regions explicitly (in the worst-case) in time $O(n^{d})$ (which is our target runtime complexity). We avoid the need to compute the Tukey regions explicitly by using an approximation scheme for estimating the volume of each region and for sampling a point from such a region, a scheme that is based on the volume estimation framework of Lovász and Vempala [13] and of Cousins and Vempala [7]. The second issue is that it might be the case that all Tukey regions for large values of $k$ are degenerate, i.e., have volume 0, in which case this strategy might fail to identify a point in the convex hull of the input points.

We remark that it is easy to construct algorithms with running time polynomial in the input size $n$, when $n$ grows exponentially in $d$. (In that case one can solve the problem using a reduction to the 1-dimensional case.) In this work we aim to reduce the running time while keeping the sample complexity $n$ at most polynomial in $d$ and in $\log|G|$.

Let $G^{*}$ denote the set of all finite databases over a grid $G$, and let $F$ be a finite set. Given a database $S\in G^{*}$, a quality (or scoring) function $q:G^{*}\times F\rightarrow{\mathbb{N}}$ assigns a number $q(S,f)$ to each element $(S,f)\in G^{*}\times F$, identified as the “quality” of $f$ with respect to $S$. We say that the function $q$ has sensitivity $\Delta$ if for all neighboring databases $S$ and $S^{\prime}$ and for all $f\in F$ we have $|q(S,f)-q(S^{\prime},f)|\leq\Delta$.

The exponential mechanism of McSherry and Talwar [15] privately identifies an element $f\in F$ with large quality $q(S,f)$. Specifically, it chooses an element $f\in F$ randomly, with probability proportional to $\exp\left({\varepsilon}\cdot q(S,f)/(2\Delta)\right)$. The privacy and utility of the mechanism are:

The Tukey depth [17] ${\sf td}_{S}(q)$ of a point $q$ with respect to $S$ is the minimum number of points of $S$ we need to remove to make $q$ lie outside the convex hull of the remaining subset. Equivalently, ${\sf td}_{S}(q)$ is the smallest number of points of $S$ that lie in a closed halfspace containing $q$. We will write ${\sf td}(q)$ for ${\sf td}_{S}(q)$ when the set $S$ is clear from the context. It easily follows from Helly’s theorem that there is always a point of Tukey depth at least $n/(d+1)$ (see, e.g., [14]). We denote the largest Tukey depth of a point by ${\sf td}_{\rm max}(S)$ (the maximum Tukey depth is always at most $n/2$).

We define the regions $D_{\geq k}(S)=\left\{q\in[0,1]^{d}\mid{\sf td}_{S}(q\geq k\right\}$ and $D_{k}(S)=D_{\geq k}(S)\setminus D_{\geq k+1}(S)$ for $k=0,\ldots,{\sf td}_{\rm max}(S)$. Note that $D_{\geq 1}$ is the convex hull of $S$ and that $D_{\geq 0}=[0,1]^{d}$. It is easy to show that $D_{\geq k}$ is convex; $D_{\geq k}$ is in fact the intersection of all (closed) halfspaces containing at least $n-k+1$ points of $S$; see [16]. It is easy to show that all this is true also when $S$ is degenerate. See Figure 1 for an illustration. The following lemma is easy to establish.

For every affine subspace $f$, of dimension $0\leq j\leq d-1$, spanned by some subset of (at least) $j+1$ points of $G$, we denote by $c(f)$ the number of points of $S$ that $f$ contains, and refer to it as the size of $f$.

We start by computing $c(f)$ for every subspace $f$ spanned by points of $S$, as follows. We construct the (potentially degenerate) arrangement ${\cal A}(S^{*})$ of the set $S^{*}$ of the hyperplanes dual to the points of $S$. During this construction, we also compute the multiplicity of each flat in this arrangement, namely, the number of dual hyperplanes that contain it. Each intersection flat of the hyperplanes is dual to an affine subspace $f$ defined by the corresponding subset of the primal points of $S$ (that it contains), and $c(f)$ is the number of dual hyperplanes containing the flat. In other words, as a standard byproduct of the construction of ${\cal A}(S^{*})$, we can compute the sizes of all the affine subspaces that are spanned by points of $S$, in $O(n^{d})$ overall time.

We define

and compute $M^{\prime}_{i}=M_{i}+Y_{i}$, where $Y_{i}$ is a random variable drawn (independently for each $i$) from a Laplace distribution with parameter $b:=\frac{1}{{\varepsilon}}$ centered at the origin. (That is, the probability density function of $Y_{i}$ is $\frac{1}{2b}e^{-|x|/b}=\frac{{\varepsilon}}{2}e^{-{\varepsilon}|x|}$.)

Our algorithm now uses a given confidence parameter $\beta\in(0,1)$ and proceeds as follows. If for every $j=0,\ldots,d-1$

we apply the base case. Otherwise, set $j$ to be the smallest index such that

Having fixed $j$, we find (privately) a subspace $f$ of dimension $j$ that contains a large number of points of $S$. To do so, let $Z_{j}$ be the collection of all $j$-dimensional subspaces that are spanned by $j+1$ affinely independent points of the grid $G$ (not necessarily of $S$). We apply the exponential mechanism to pick an element of $Z_{j}$, by assigning a score $s(f)$ to each subspace of $Z_{j}$, which we set to be

if $j\geq 1$, and $s(f)=c(f)$ if $j=0$. Note that by its definition, $s(f)$ is zero if $f$ is not spanned by points of $S$. Indeed, in this case the $c(f)$ points contained in $f$ span some subspace of dimension $\ell\leq j-1$ and therefore $M_{j-1}$ must be at least as large as $c(f)$. We will shortly argue that $s(f)$ has sensitivity at most $2$ (Lemma 4.4), and thus conclude that this step preserves the differential privacy of the procedure.

We would like to apply the exponential mechanism as stated above in time proportional to the number of subspaces of non-zero score, because this number depends only on $n$ (albeit being exponential in $d$) and not on (the much larger) $X$. However, to normalize the scores to probabilities, we need to know the number of elements of $Z_{j}$ with zero score, or alternatively to obtain the total number of subspaces spanned by $j+1$ points of $G$ (that is, the size of $Z_{j}$).

We do not have a simple expression for $|Z_{j}|$ (although this is a quantity that can be computed, for each $j$, independently of $S$, once and for all), but clearly $|Z_{j}|\leq X^{d{(j+1)}}$. We thus augment $Z_{j}$ (only for the purpose of analysis) with $X^{d{(j+1)}}-|Z_{j}|$ “dummy” subspaces, and denote the augmented set by $Z_{j}^{\prime}$, whose cardinality is now exactly $X^{d{(j+1)}}$. We draw a subspace $f$ from $Z_{j}^{\prime}$ using the exponential mechanism. To do so we need to compute, for each score $s\geq 0$, the number $N_{s}$ of elements of $Z^{\prime}_{j}$ that have score $s$, give each such element weight $e^{{\varepsilon}s/4}$, choose the index $s$ with probability proportional to $N_{s}e^{{\varepsilon}s/4}$, and then choose uniformly a subspace from those with score $s$. It is easy to check that this is indeed an implementation of the exponential mechanism as described in Section 2.1.

If the drawing procedure decides to pick a subspace that is not spanned by points of $S$, or more precisely decides to pick a subspace of score $0$, we stop the whole procedure, with a failure. If the selected score is positive, the subspace to be picked is spanned by $j+1$ points of $S$, and those subspaces are available (from the dual arrangement construction outlined above). We thus obtain a selected subspace $f$ (by randomly choosing an element from the available list of these subspaces), and apply the algorithm recursively within $f$, restricting the input to $S\cap f$. (Strictly speaking, as noted above, we apply the algorithm to a projection of $f$ onto a suitable $j$-dimensional coordinate subspace.)

It follows that we can implement the exponential mechanism on all subspaces in $Z_{j}^{\prime}$ in time proportional to the number of subspaces spanned by points of $S$, which is $O(n^{d})$, and therefore the running time of this subspace selection procedure (in each recursive call) is $O(n^{d})$.