# Instance-Optimal Private Density Estimation in the Wasserstein Distance

Vitaly Feldman
Apple  Audra McMillan
Apple  Satchit Sivakumar111Work partially done while author was an intern at
Apple.
Boston University  Kunal Talwar
Apple

###### Abstract

Estimating the density of a distribution from samples is a fundamental problem
in statistics. In many practical settings, the Wasserstein distance is an
appropriate error metric for density estimation. For example, when estimating
population densities in a geographic region, a small Wasserstein distance
means that the estimate is able to capture roughly where the population mass
is. In this work we study differentially private density estimation in the
Wasserstein distance. We design and analyze instance-optimal algorithms for
this problem that can adapt to easy instances.

For distributions $P$ over $\mathbb{R}$, we consider a strong notion of
instance-optimality: an algorithm that uniformly achieves the instance-optimal
estimation rate is competitive with an algorithm that is told that the
distribution is either $P$ or $Q_{P}$ for some distribution $Q_{P}$ whose
probability density function (pdf) is within a factor of 2 of the pdf of $P$.
For distributions over $\mathbb{R}^{2}$, we use a different notion of instance
optimality. We say that an algorithm is instance-optimal if it is competitive
with an algorithm that is given a constant-factor multiplicative approximation
of the density of the distribution. We characterize the instance-optimal
estimation rates in both these settings and show that they are uniformly
achievable (up to polylogarithmic factors). Our approach for $\mathbb{R}^{2}$
extends to arbitrary metric spaces as it goes via hierarchically separated
trees. As a special case our results lead to instance-optimal private learning
in TV distance for discrete distributions.

###### Contents

  1. 1 Introduction
    1. 1.1 Our Results
    2. 1.2 Techniques
      1. 1.2.1 Distributions over $\mathbb{R}$:
      2. 1.2.2 Distributions on HSTs
  2. 2 Preliminaries
    1. 2.1 Differential Privacy
  3. 3 On Instance Optimality
    1. 3.1 Local Estimation Rates
    2. 3.2 Locally Minimal Algorithms
    3. 3.3 Relaxed Definitions
  4. 4 Additional Related Work
  5. 5 Distribution Estimation on Hierarchically Separated Trees
    1. 5.1 Preliminaries on Hierarchically Separated Trees
    2. 5.2 The Target Estimation Rate
      1. 5.2.1 Reduction to Estimation in TV distance of Discrete Distributions
      2. 5.2.2 Characterizing Target Estimation Rate for Discrete Distributions
    3. 5.3 An $\varepsilon$-DP Distribution Estimation Algorithm
  6. 6 Instance Optimal Density Estimation on $\mathbb{R}$ in Wasserstein distance
    1. 6.1 General Lower Bound
      1. 6.1.1 The Privacy Term
      2. 6.1.2 Empirical Term
    2. 6.2 Upper Bound
      1. 6.2.1 Algorithm for density estimation
  7. A Preliminaries
    1. A.1 Distribution Distances
    2. A.2 Differential Privacy
  8. B Experiment Details
  9. C Appendix for Section 5
  10. D Local Minimality in the High Dimensional Setting
  11. E Differentially Private Quantiles
  12. F Proofs in Section 6
    1. F.1 Omitted Proofs in Section 6.1.1
    2. F.2 Omitted proofs in Section 6.1.2
    3. F.3 Omitted Proofs in Section 6.2
    4. F.4 Local Minimality in the One-Dimensional Setting

##  1 Introduction

Distribution estimation is a fundamental problem in statistics. In this work,
we focus on the problem of learning the density of a distribution over a low-
dimensional real space. Our motivation for studying this problem comes from
practical problems such as estimating the population density in a geographical
area (defined by bounded two dimensional space, for e.g. $[0,\ell]^{2}$),
learning the distribution of accuracy of a machine learning model (i.e. a
distribution over $[0,1]$), estimating the average temperature across
latitude, longitude, and altitude (i.e. a distribution over $[0,\ell]^{3}$)
etc.

In this work, we are interested in the non-parametric version of this
question, where we make no assumptions on the form of the distribution we are
learning. This is frequently of interest in practice, where population
densities for example may change over time (become more or less concentrated),
and it is difficult to specify a meaningful parametric class that will
simultaneously capture all densities of interest. Given estimation is often
done using sensitive data (for e.g. health data), our interest in this
question is in, and consequently all our results are for, the differentially
private version of this question. While we believe our results in the non-
private setting are also novel and interesting, we view the private results as
our main contribution.

Any statistical algorithm learning from samples is inexact. The appropriate
gauge to measure the (in)accuracy of a density estimation algorithm depends on
how this density estimate is used. In this work, we focus on the Wasserstein
distance between the original distribution and the learnt distribution as our
measure of accuracy. Known by many names (Earthmover distance, Kantorovich
distance, Optimal Transport distance), this distance is defined over any
distance metric $d$ as the minimum over all couplings $\pi$ from $P$ to $Q$ of
the quantity $\operatorname*{\mathbb{E}}_{x\sim P}[d(x,\pi(x))]$. It is
arguably one of the most natural ways to define distances between
distributions over a metric space and has been extensively studied (see
Section 4) . We note that Wasserstein distance is particularly salient in many
practical applications of density estimation where the geometry of the space
is significant. As a simple example, when creating population density
estimates, if the population is concentrated in a few cities, then outputting
a distribution concentrated close to these cities (even if not exactly at the
cities) is intuitively better than outputting a distribution that is more
spread out. Metrics such as TV distance that do not incorporate the geometry
of the space do not capture this nuance. Additionally, Wasserstein distance is
versatile and can be adapted to the setting of interest by varying the metric.
In the case of the metric being a discrete metric with
$d(x,y)=\mathbf{1}(x\neq y)$, it reduces to the commonly used total variation
$(\text{\rm TV})$ distance. Our focus in this work is on the case of Euclidean
distance metric on $[0,1]$ or $[0,1]^{2}$, though our results apply to both to
higher-dimensional Euclidean space as well as to any finite metric. In the
$[0,1]$ case (with the standard Euclidean metric), the Wasserstein distance is
equivalent to the total area between the cumulative distribution functions.

The problem of learning a distribution under Wasserstein distance has a long
history, starting with [Dud69] proving worst-case bounds on the rate of
convergence of the Wasserstein distance between the empirical distribution
$\hat{P}_{n}$ and the target distribution $P$ over $\mathbb{R}^{d}$.
Similarly, this question for the case of the discrete metric
($d(x,y)=\mathbf{1}(x\neq y)$) has been very well studied. However, most known
results for this problem look at it from the point of view of worst-case
analysis. This can paint a rather pessimistic picture. For example, the
minimax rate of $\varepsilon$-privately learning a discrete distribution over
$\\{0,\dots,k\\}$ in TV distance (i.e. Wasserstein with the discrete metric
described above) scales linearly with $k$, which can be prohibitive for large
support size $k$. For Wasserstein distance with $\ell_{2}$ norm, the rate of
convergence of the empirical distribution suffers a curse of dimensionality,
with the worst-case error between the distribution and the empirical
distribution being $\Theta(n^{-\frac{1}{d}})$ for distributions over
$[0,1]^{d}$. For the differentially private version of this question, recent
works [BSV22, HVZ23] have shown that the optimal Wasserstein minimax error
between the sample and the private estimate is $\tilde{\Theta}((\varepsilon
n)^{-\frac{1}{d}})$. This worst-case analysis viewpoint fails to distinguish
between algorithms that perform very differently on the types of instances one
may see in practice. In particular, many practical distributions may be more
feasible to estimate than suggested by the minimax rate. As an example, Figure
1 shows the cumulative distribution function of a bimodal distribution on
$[0,1]$ with very sparse support, and the cdf learnt by a minimax optimal
algorithm, as well as an algorithm we present in this work . As is clear from
the figure, the minimax optimal algorithm is easily outperformed. This
phenomenon only gets worse in higher dimensions. Similarly, if the
distribution in $\real^{d}$ lies on a $k$-dimensional subspace, the worst-case
error scaling with $\tilde{O}((\varepsilon n)^{-\frac{1}{d}})$ is
significantly larger than our algorithm’s scaling of $\tilde{O}((\varepsilon
n)^{-\frac{1}{k}})$.

Figure 1: (Left) A sparsely supported distribution on integers [0,999] (pdf).
(Right) CDF for the same distribution (green, solid line), along with a (non-
private) minimax optimal learnt distribution (blue, dashed line), as well as
1-DP instance-optimal algorithm (red, dotted), both learnt from the same 1600
samples. The $W_{1}$ error for the minimax optimal algorithm is 13.4, whereas
the DP estimated distribution has $W_{1}$ error of $0.86$. While this example
is artificial, it demonstrates the large potential gap between minimax optimal
and instance optimal algorithms on specific instances.

This motivates the problem of viewing this question through the lens of
instance optimality. 222c.f. related work section for discussion of other
beyond worst case analysis approaches for this question and Section 3 for a
more in-depth discussion of our approach. Briefly, instance optimal algorithms
are those that on any given instance of the problem, are able to perform
competitively with what any algorithm can do on this instance. Let
$\mathcal{M}$ be a class of algorithms of interest (e.g. all
$(\varepsilon,\delta)$-differentially private algorithms) and $cost(\cdot,P)$
be a cost measure for an instance $P$. In our setting, we have a distribution
$P$ over a metric space, and given a set $\hat{P}_{n}$ of $n$ samples from
$P$, we want to learn an estimate $\mathcal{A}(\hat{P}_{n})$ for the
distribution. Our measure of performance is the Wasserstein distance
$\mathcal{W}$, so
$cost(\mathcal{A},P)=\mathbb{E}[\mathcal{W}(P,\mathcal{A}(\hat{P}_{n}))]$. We
would ideally like to say that an algorithm $\mathcal{A}$ is $\alpha$-instance
optimal in a class $\mathcal{M}$ if for all instances $P$, and all
$\mathcal{A^{\prime}}\in M$,

$\displaystyle cost(\mathcal{A}(\hat{P}_{n})),P)\leq\alpha\cdot
cost(\mathcal{A^{\prime}}(\hat{P}_{n})),P).$ (InstanceOptimality-Ideal)

The reader would have noticed that this definition is however impossible to
achieve except for trivial classes $\mathcal{M}$. The algorithm
$\mathcal{A}^{\prime}$ that ignores its input and always outputs $P$ makes the
right hand side 0. However, this algorithm performs poorly on any
distributions far from $P$ and so is not a reasonable benchmark. A common
approach in many works is to measure the performance of the competing
algorithm $\mathcal{A}^{\prime}$ not just on the given instance, but on a
small neighborhood around it. Thus we say that that an algorithm $\mathcal{A}$
is $\alpha$-instance optimal amongst a class $\mathcal{M}$ with respect to a
neighborhood function $\mathcal{N}$ if for all instances $P$, and all
$\mathcal{A^{\prime}}\in M$

$\displaystyle
cost(\mathcal{A}(\hat{P}_{n})),P)\leq\alpha\cdot\sup_{P^{\prime}\in\mathcal{N}(P)}cost(\mathcal{A^{\prime}}(\hat{P^{\prime}}_{n})),P^{\prime}).$

In other words, the benchmark we evaluate against is the cost of the best
algorithm for a neighborhood $\mathcal{N}(P)$ _that knows this neighborhood_.
We would like our algorithm $\mathcal{A}$, that is not tailor-made for
$\mathcal{N}(P)$, to nevertheless be competitive against this benchmark.

This definition is general, and captures most notions of instance optimality
that have been studied in the literature. The set $\mathcal{N}(P)$ must be
carefully defined for this notion to be meaningful; we can always define
$\mathcal{N}(P)$ to be the set of all instances whence this notion reduces to
worst-case analysis. In many previous works, this neighborhood map has been
defined to capture the belief that any natural algorithm must not have
significantly different performance on different members of $\mathcal{N}(P)$.
For example, [FLN01, ABC17, VV16, OS15, GKN20] include in $\mathcal{N}(P)$
appropriate renamings of $P$ to capture some kind of permutation invariance of
natural algorithms. In statistics, one often enforces that the cardinality of
$\mathcal{N}(P)$ is 2, often called the _hardest one-dimensional subproblem_
[CL15, AD20, DLSV23]. Some recent works in privacy [HLY21, DKSS23] have
defined instance optimality w.r.t. neighboring datasets obtained by deleting a
small number of data points. Any reasonable definition of instance optimality
for a problem must justify its choice of the neighborhood map; similar choices
must be justifiable in every other notion of beyond worst case analysis
[Rou21]. In instance-optimality definitions, this choice of neighborhood is
what encapsulates what class of domain-specific algorithms our algorithm
competes against. A good definition thus depends on the context and on the
kind of domain knowledge we imagine an expert designing a custom algorithm for
an application may have. Ideally, the definition is broad (i.e. the
neighborhoods $\mathcal{N}(P)$ are sufficiently contained) so that in a large
class of applications, we expect the domain knowledge to not be enough to rule
out any member of $\mathcal{N}(P)$. We discuss this general definition of
instance optimality further in Section 3. We remark that for reasonable
neighborhood maps, this is an extremely strong requirement: an instance-
optimal algorithm must simultaneously do well on every single input, in fact
as well as any other algorithm that is given this neighborhood
$\mathcal{N}(P)$ in advance!

Instance optimality guarantees are most useful when there is a big difference
between achievable utility guarantees for typical cases and the worst-case
utility guarantees. Wasserstein estimation is an example of such a problem. We
will see that achievable utility bounds for, for example, concentrated
distributions are a lot better than worse case distributions. Our definition
of instance optimality is particularly suitable for metric spaces, and our
notion of neighborhood allows the target utility bound to adapt to the
distribution. We note that for estimation in Wasserstein distance with
practically important metrics such as $\ell_{1}$ and $\ell_{2}$ norms, it is
unclear if existing instance optimality definitions (using notions of
neighborhood discussed above) capture this. For example, for discrete
distributions, setting the neighborhood to be all permutations of the
distribution destroys all structure of the distribution (for e.g.
concentration), and hence performance on this neighborhood may not capture the
relative ease of estimation of a concentrated distribution. Similar problems
apply to other previously studied definitions of instance optimality, which
are not well-suited to density estimation with error metrics that incorporate
the geometry of the metric space. See Section 3 and Section 4 for further
discussion on the inadequacy of existing instance optimality definitions for
our setting of interest.

Our notion of neighborhood will correspond to small balls in one of the
strictest notions of distance between distributions. Recall that for
distributions $P,Q$ on $X$, $D_{\infty}(P,Q)$ is defined as $\sup_{x\in
X}\max\left(\ln\frac{P(x)}{Q(x)},\ln\frac{Q(x)}{P(x)}\right)$. Our
neighborhood map $\mathcal{N}$ will have the property that for all $P$, and
for all $Q\in\mathcal{N}(P)$, $D_{\infty}(P,Q)\leq\ln 2$. This corresponds to
the benchmark algorithm $\mathcal{A}^{\prime}$ being given as auxiliary input
a multiplicative constant factor approximation to the probability density
function $P(x)$ (and we can replace the constant $2$ by any constant). In
particular, an algorithm that knows the support of the distribution $P$ will
not be able to do much better than our algorithm that gets no such
information. Notice that this implicitly implies that our algorithm is able to
exploit sparsity in the data distribution since it is competitive with an
algorithm that is told the support. In the one-dimensional real case we can
achieve an even stronger notion of instance-optimality. In this case
$\mathcal{N}(P)$ is defined to be $\\{P,Q\\}$ where $Q$ is a distribution with
$D_{\infty}(P,Q)\leq\ln 2$. This is a strengthening of the rate defined by the
hardest one-dimensional subproblem.

We also give a definition that captures another aspect of instance optimality,
related to the notion of super efficiency, that we term local minimality in
Section 3. Informally, local minimality says that if any comparator algorithm
does better than $\mathcal{A}$ on $P$, then there is a distribution $Q$ in the
neighborhood of $P$ where $\mathcal{A}$ does better than the comparator.
Approximate local minimality relaxes the latter condition to being better than
some constant times the comparator. The two definitions of approximate local
minimality and instance optimality are in general incomparable (see Section 3)
but for suitable smooth algorithms, we show that these definitions are
equivalent. Our algorithms, both for the 1-dimensional and the case of general
metric spaces approximately satisfy both these definitions.

In order to show that the instance optimality definition is achievable, we
give both algorithmic upper bounds and matching, up to logarithmic factors,
theoretical lower bounds. The algorithms we use in our upper bounds are built
largely from ingredients previously used for similar problems. We see this as
an asset since these algorithms are implementable in practice. A key
ingredient that we do introduce is the use of randomised HST approximation of
finite metric spaces. This replaces deterministic hierarchical decompositions
that were used in prior work, allowing us to gain tighter utility guarantees.
Our main conceptual contribution is to introduce what we believe to the right
notion of instance optimality for this problem, including the definition of a
meaningful neighbourhood function. The main technical challenge is in the
lower bounds, which require carefully building nets of distributions within
each neighborhood $\mathcal{N}(P)$ that allow us to use a slight
generalisation of DP Assoud’s Lemma to give a lower bound on the target
estimation rate for each distribution $P$.

###  1.1 Our Results

Given an estimation algorithm $\mathcal{A}:\mathcal{X}^{n}\to\mathcal{M}$, the
estimation rate of $\mathcal{A}$ for distribution $P$ is:

$\displaystyle R_{\mathcal{A},n}(P)=\inf_{t\in\mathbb{R}}$
$\displaystyle\\{t:\text{w.p. $\geq 0.75$ over ${\bf x}\sim P^{n}$ and the
randomness of the algorithm, }\mathcal{W}_{d}(\mathcal{A}({\bf x}),P)\leq
t\\}.$ (1)

We start by stating an informal version of our result in the one-dimensional
real case.

######  Theorem 1.1 (Informal 1-dimensional result).

Let $\varepsilon,\gamma\in(0,1]$. There is an $\varepsilon$-differentially
private algorithm $\mathcal{A}$ such that, for all distributions $P$ supported
in $[0,1]$, for all natural numbers
$n>\frac{\operatorname{polylog}1/\gamma}{\varepsilon}$, there exists a
distribution $Q$ (with $D_{\infty}(P,Q)\leq\ln 2$) such that the following is
satisfied.

For any $\varepsilon$-DP algorithm $\mathcal{A}^{\prime}$, with probability at
least $0.75$ over the randomness of ${\bf x}\sim P^{n}$ and additional
randomness of the algorithm,

$\displaystyle\mathcal{W}(\mathcal{A}(\hat{P}_{n}),P)$
$\displaystyle\leq\operatorname{polylog}n\cdot\sup_{P^{\prime}\in\\{P,Q\\}}R_{\mathcal{A}^{\prime},n^{\prime}}(P^{\prime})+\gamma$

where $n^{\prime}\approx\frac{n}{\operatorname{polylog}n/\gamma}$

In this one-dimensional case, our algorithm is based on DP quantile
estimation. The additive $\gamma$ term can be made polynomially small. The
lower bound is based on (differentially private) simple hypothesis testing
where for each distribution $P$, we find a distribution in $\mathcal{N}(P)$
that is indistinguishable from $P$ given $n$ samples but also sufficiently far
from $P$ in Wasserstein distance.

Extending the quantiles based approach from the one dimensional setting to
even the two dimensional setting is challenging, as there is no “right” way to
generalize quantiles to dimensions 2 or beyond. Several previous works on
Wasserstein density estimation (e.g. [BNNR09]) have used a hierarchical
decomposition approach to address this question. A hierarchical approach has
also been used in various more practical works on private density estimation
(e.g. [CB22, QYL12, BKM+21, MJT+22, ZXX16]). These works focus on practical
performance and do not offer tight theoretical bounds. A hierarchical approach
was also used by [GHK+23], who proved theoretical bounds for a related
problem, but not through the lens of instance optimality. We compare our
results to theirs in more detail later in this section.

The use of deterministic hierarchical decompositions in all these papers means
that some points that are very close (but on opposite sides of the boundaries
of the hierarchical decomposition) get mapped to relatively far points,
resulting in high distortion factors that are not appropriate for instance
optimality.

Inspired by the above approaches but noting their constraints, we use a
randomized embedding into hierarchically separated trees instead of a
deterministic one. We define our algorithm on any hierarchically separated
tree metric and use the fact that there is a randomized embedding of
$[0,1]^{2}$ on a hierarchically separated tree metric space with low
distortion. This, along with some other important technical modifications
(such as truncating low values to $0$), allows us to analyze a variant of the
above practical algorithms theoretically and show that it satisfies our strong
notion of instance optimality, up to polylogarithmic factors in the number of
samples.

######  Theorem 1.2 (Informal two-dimensional result).

There is a polynomial time $\varepsilon$-differentially private algorithm
$\mathcal{A}$ that for any distribution $P$ on $[0,1]^{2}$, any integer $n$,
and any $\varepsilon$-DP algorithm $\mathcal{A}^{\prime}$ with probability at
least 0.75, satisfies

$\displaystyle\mathcal{W}_{2}(\mathcal{A}(\hat{P}_{n}),P)$
$\displaystyle\leq(\log
n)^{O(1)}\sup_{P^{\prime}:D_{\infty}(P,P^{\prime})\leq\ln
2}\operatorname*{\mathbb{E}}[\mathcal{W}_{2}(\mathcal{A}^{\prime}(\hat{P^{\prime}}_{n^{\prime}}),P^{\prime})],$

where $n^{\prime}\approx\frac{n}{\operatorname{polylog}n}$. Here, the
expectation is taken over the internal coin tosses of $\mathcal{A}$ as well as
over the choice of the i.i.d. samples $\hat{P}_{n}$.

In fact, since our algorithm is defined on any hierarchically separated tree
metric space, it has the added bonus of giving instance optimality results for
any finite metric space (since powerful results [Bar96, FRT03] show that any
finite metric space can be embedded in a hierarchically separated tree metric
space with a distortion factor at most logarithmic in the size of the metric
space).

######  Theorem 1.3 (Informal finite metric result).

Let $(\mathcal{X},d)$ be an arbitrary metric space with diameter $1$. There is
a polynomial time $\varepsilon$-differentially private algorithm $\mathcal{A}$
such that for any distribution $P$ on $X$ any integer $n$ and any
$\varepsilon$-DP algorithm $\mathcal{A}^{\prime}$ with probability at least
0.75, satisfies

$\displaystyle\mathcal{W}(\mathcal{A}(\hat{P}_{n}),P)$
$\displaystyle\leq(\log|\mathcal{X}|\cdot\log
n)^{O(1)}\sup_{P^{\prime}:D_{\infty}(P,P^{\prime})\leq\ln
2}\operatorname*{\mathbb{E}}[\mathcal{W}(\mathcal{A}^{\prime}(\hat{P^{\prime}}_{n^{\prime}}),P^{\prime})],$

where $n^{\prime}\approx\frac{n}{\operatorname{polylog}n}$. Here, the
expectation is taken over the internal coin tosses of $\mathcal{A}$ as well as
over the choice of the i.i.d. samples $\hat{P}_{n}$.

Our lower bound result is actually slightly stronger than stated in Theorem
1.3 since it holds not only for $\varepsilon$-DP, but also for
$(\varepsilon,\delta)$-DP. At this point, we also compare specifically to the
paper of [GHK+23] who give an algorithm for obtaining two-dimensional heatmaps
and analyze it theoretically. They focus on the empirical version of a variant
of this problem as opposed to the population version, and aim to compete with
the best $k$-sparse distribution. Their algorithm takes the sparsity parameter
$k$ as input in order to set parameters and achieves additive error
$\sqrt{k}/n$ (and a constant multiplicative factor). On the other hand, our
algorithm also performs better for sparse distributions but is automatically
adaptive to the sparsity (and hence doesn’t need to take it as an input).
Additionally the additive term in our work can be made polynomially small (for
any polynomial) in $n$ at a logarithmic cost to the multiplicative error
(regardless of the sparsity of the distribution). On the other hand, for large
$k$ their results have additive error that scales with $1/\sqrt{n}$. Their use
of a deterministic hierarchical decomposition makes their algorithm unsuitable
for our notion of instance optimality (as discussed earlier), and it is
unclear if their algorithm can be directly extended to all finite metric
spaces.

Note that instance optimality for all finite metric spaces implies instance
optimality results for a wide variety of applications not addressed in prior
work. For example, our results immediately extend to other low-dimensional
real spaces with arbitrary metrics (for e.g. $\ell_{p}$ norms). They also give
non-trivial improvements on worst-case analysis for higher-dimensional spaces
that are not the main focus of our work (for $[0,1]^{d}$, we can use a fine
grid of size $(1/(\eta/\sqrt{d}))^{d}$ at an additive cost of $\eta$ in the
Wasserstein distance in order to create a finite metric space to apply our
result on. Since the dependence on $|\mathcal{X}|$ in the result above is
logarithmic, this translates to a $d\log\frac{d}{\eta}$ multiplicative
overhead term replacing the $\log|\mathcal{X}|$ factor above. While this is
still a significant overhead, all previous results on density estimation in
the Wasserstein distance (in both the private and non-private literature) are
worst case, where the sample complexity is exponential in $d$. Since our
results only have a polynomial dependence in $d$ over the optimal error, this
is a non-trivial improvement over worst-case error, even when $d$ is large.

Another immediate application of our results is to give (to the best of our
knowledge) new bounds for private estimation of discrete distributions in TV
distance. Generally, for learning a discrete distribution defined by
probabilities $\\{p_{1},\ldots,p_{k}\\}$, our results lead to a rate (up to
polylogarithmic factors) of
$\sum_{i}\min\Big{\\{}p_{i}(1-p_{i}),\sqrt{\frac{p_{i}(1-p_{i})}{n}}\Big{\\}}+\sum_{i}\min\left(p_{i},(1-p_{i}),\frac{1}{\varepsilon
n}\right).$

This can give significant improvements over the worst case bounds for
practically important distributions. The minimax rate is linear in the support
size $k$, namely $\Theta(k/\varepsilon n)$ (for sufficiently small
$\varepsilon$). Now, consider the following power-law distribution over
support size $k$: $p(i)\propto i^{-2}$. (Power law distributions arise
frequently in practice for e.g. frequencies of family names, sizes of power
outages etc. all follow power law distributions.) Applying our result above
gives a bound that is $\tilde{O}\left(\min\\{\frac{k}{\varepsilon
n},\frac{1}{\sqrt{\varepsilon n}}\\}\right)$, which is much better than the
worst case bound for large support distributions.

Our result also applies to other practically important settings such as
building lists of popular sequences such as n-grams over words. We leave open
the questions of designing instance-optimal algorithms for other practically
important questions in private learning and statistics, and of designing
better instance optimal algorithms for higher dimensional spaces. We also
leave open the question of removing the polylogarithmic factors in our
instance optimality bounds.

###  1.2 Techniques

####  1.2.1 Distributions over $\mathbb{R}$:

We start by describing the rate we obtain for distributions $P$ over
$\mathbb{R}$.In order to state the rate, we will use $q_{\alpha}$ to represent
the $\alpha$-quantile of the distribution $P$ and use $P|_{a,b}$ to define a
certain restricted distribution described below. The rate consists of three
terms and roughly looks as follows— we suppress logarithmic factors in $n$.

$R_{\mathcal{A},n}(P)=\tilde{O}\left(\mathbb{E}\left[\mathcal{W}\left(P,\hat{P}_{n}\right)\right]+\frac{1}{\varepsilon
n}\left(q_{1-\frac{1}{\varepsilon n}}-q_{\frac{1}{\varepsilon
n}}\right)+\mathcal{W}(P,P|_{q_{\frac{1}{\varepsilon
n}},q_{1-\frac{1}{\varepsilon n}}})\right),$

The first term is $\mathbb{E}[\mathcal{W}(P,\hat{P}_{n})]$, the expected
Wasserstein distance between the true distribution and the empirical
distribution over $n$ samples, and is the non-private term. The remaining two
terms represent the cost of privacy- the first is a specific interquantile
distance, roughly $\frac{1}{\varepsilon n}(q_{1-\frac{1}{\varepsilon
n}}-q_{\frac{1}{\varepsilon n}})$, and the second can be thought of as
capturing the weight of the tails- represented by the Wasserstein distance
between $P$ and a ‘restricted’ version of $P$ with its tails chopped off (i.e.
the cumulative distribution function is $0$ below $q_{1/\varepsilon n}$ and
$1$ above $q_{1-1/\varepsilon n}$ and identical to $P$ otherwise). Observe
that all $3$ of the terms above are smaller for distributions with small
support or greater concentration, and hence the rate adapts to the hardness of
the distributions.

Upper Bounds: The upper bound involves estimating roughly $\varepsilon n$
equally spaced quantiles of the empirical distribution differentially
privately (using a known private CDF estimation algorithm), and placing
roughly $1/\varepsilon n$ mass at each of the estimated quantile points. For
the analysis, the intuition for each of the terms is as follows: since we only
have access to the empirical distribution, the non-private term
$\mathbb{E}[\mathcal{W}(P,\hat{P}_{n})]$ comes from that. Next, if the
quantile estimates are good, then the pointwise CDF differences between the
empirical distribution and the estimated distribution are at most
$1/\varepsilon n$ (due to the discretization), and so we will pay
$1/\varepsilon n$ multiplied by the interquantile distance of the empirical
distribution. This aligns with the accuracy of state-of-the-art DP quantile
estimation algorithms. Finally, since the distribution is restricted to the
estimated quantiles, the distribution is $0$ before the first estimated
quantile and $1$ above the last estimated quantile and so we pay the
Wasserstein distance between the empirical distribution and a restricted
version of the empirical distribution. Some care needs to be taken while
reasoning about expectation versus high probability (for various terms), and
in relating population quantities to empirical quantities (which we do using
various concentration inequalities). Details can be found in Section 6.2.

Lower Bounds: We prove that the private and non-private terms are lower bounds
separately. Both proofs follow the same framework. The idea is that given
knowledge of two distributions $P$ and $Q$, we can use a (private) Wasserstein
estimation algorithm to construct a hypothesis test distinguishing $P$ from
$Q$. If the (private) estimate for $P$ and $Q$ with $n$ samples gives error
smaller than $\frac{1}{2}\mathcal{W}(P,Q)$, we can use this to distinguish $P$
from $Q$. This would give a contradiction if $P$ and $Q$ are (privately)
indistinguishable with $n$ samples. Hence, this would give a lower bound of
$\frac{1}{2}\mathcal{W}(P,Q)$ on the error of the Wasserstein estimation
algorithm on $P$ or $Q$.

Thus the task reduces to constructing a distribution $Q$ that satisfies three
properties: 1) it is (privately) indistinguishable from $P$ given $n$ samples,
2) the Wasserstein distance between $P$ and $Q$ is sufficiently large, 3)
$D_{\infty}(P,Q)\leq\ln 2$. The main technical work is in identifying a
distribution $Q$ that satisfies these properties.

For the privacy term, we construct the distribution $Q$ by taking half the
mass from the first $1/\varepsilon n$-quantile of $P$ (scaling the density
function by half) and moving it to the last $1/\varepsilon n$-quantile of $P$
(scaling the density function by $3/2$). The third property is satisfied by
definition, so we reason about the other two. Intuitively, since the
Wasserstein distance captures how hard it is to ‘move’ $P$ to $Q$, this mass
needs to move at least the interquantile distance to change $P$ to $Q$. This
implies that the Wasserstein distance is at least the interquantile distance
scaled by $1/\varepsilon n$, as described in the rate. Additionally, mass that
is further out in the tail needs to move more, which is captured by the
Wasserstein distance between the distribution $P$ and its ‘restriction’.
Hence, the Wasserstein distance between $P$ and $Q$ is lower bounded by these
two terms of interest. The intuition behind Property 2 is that it is hard for
any $\varepsilon$-DP algorithm to pinpoint the location of an
$\frac{1}{\varepsilon n}$-fraction of the points in the dataset. Overall, this
shows the privacy lower bound.

The non-private lower bound requires a more careful construction of $Q$. We
divide $P$ into various scales and carefully adjust them differently to obtain
the desired properties. Formally, to construct $Q$ from $P$, we consider
$q_{1/2}$ and all quantiles of the form $q_{1/2^{i}}$ and $q_{1-1/2^{i}}$ for
$i>1$. For $1\leq i<\log n$, we add mass to $[q_{1/2^{i+1}},q_{1/2^{i}})$, by
setting the density $f_{Q}$ to be $(1+\sqrt{2^{i}/n})f_{P}$ and balance out
the extra mass by setting $f_{Q}$ to be $(1-\sqrt{2^{i}/n})f_{P}$ between
$[q_{1-1/2^{i}},q_{1-1/2^{i+1}})$. For $i\geq\log n$ (i.e. the tail), we add
mass to $[q_{1/2^{i+1}},q_{1/2^{i}})$, by setting $f_{Q}$ to be
$(1+\frac{1}{2})f_{P}$ and balance out the extra mass by setting $f_{Q}$ to be
$(1-\frac{1}{2})f_{P}$ between $[q_{1-1/2^{i}},q_{1-1/2^{i+1}})$.

The third property is again trivially satisfied. For the first property,
observe that to ‘move’ $P$ to $Q$ the extra $\frac{1}{\sqrt{2^{i}n}}$ mass
between $[q_{1/2^{i+1}},q_{1/2^{i}})$ has to ‘travel’ between $q_{1/2^{i}}$
and $q_{1-1/2^{i}}$, and so the Wasserstein distance between $P$ and $Q$ can
be lower bounded by a sum of various scaled interquantile distances. We
attempt to upper bound the expected Wasserstein distance between $P$ and
$\hat{P}_{n}$ by a similar term. It is more intuitive to reason about this
using an alternative (equivalent) formulation of Wasserstein distance as the
area between the CDF curves of $P$ and $Q$. The intuition is that the expected
pointwise CDF difference between $P$ and $\hat{P}_{n}$ in the interval
$[q_{1/2^{i+1}}$, $q_{1/2^{i}})$ would be roughly $\frac{1}{\sqrt{2^{i}n}}$
(by properties of a Binomial) and hence the contribution of this interval to
the area would be roughly
$\frac{1}{\sqrt{2^{i}n}}\left(q_{1/2^{i}}-q_{1/2^{i+1}}\right)$ and similarly
for the corresponding interval $[q_{1-1/2^{i}}$, $q_{1/2^{i+1}})$. Hence, the
expected Wasserstein distance would be a sum of these scaled quantile interval
distances. We formalize this intuition using a result of Bobkov and Ledoux
[BL19] that characterizes the expected Wasserstein distance between $P$ and
$\hat{P}_{n}$ as an integral of a function of the CDF of $P$. We now have a
bound in terms of the sum of scaled quantile interval distances, but we want
to bound it by a sum of scaled interquantile distances. We can telescope the
sum to indeed bound it by a sum of scaled interquantile distances. This
establishes that $\mathcal{W}(P,Q)\geq\mathbb{E}[\mathcal{W}(P,\hat{P}_{n}]$.
Next, we show that $P$ is indistinguishable from $Q$ by analyzing the KL
divergence between $P$ and $Q$. The main idea is that high density intervals
are modified by a small multiplicative factor of roughly
$1+\frac{1}{\sqrt{n}}$, but low density intervals (with mass less than $1/n$)
are modified by a constant multiplicative factor, so overall the contribution
of each interval to the KL divergence is sufficiently small. This establishes
indistinguishability with $n$ samples. For formal details we refer the reader
to Section 6.1.

####  1.2.2 Distributions on HSTs

Since the main technical challenge of proving Theorem 1.3 is proving the
equivalent result for distributions on HST metric spaces, we focus on that
problem in this section. Standard results on low distortion embeddings of
metric spaces into HST metric spaces can be used to translate the HST result
to $[0,1]^{2}$ and to general metric spaces $X$ with $\log|X|$ overhead.

######  Definition 1.4 (Hierarchically Separated Tree).

A hierarchically separated tree (HST) is a rooted weighted tree such that the
edges between level $\ell$ and $\ell-1$ all have the same weight (denoted
$r_{\ell}$) and the weights are geometrically decreasing so
$r_{\ell+1}=(1/2)r_{\ell}$. Let $0pt$ be the depth of the tree.

HSTs can be defined with any geometric scaling but we will only need a factor
of 2 in this work. HSTs may also have arbitrary degree. A HST defines a metric
on its leaf nodes by defining the distance between any two leaf nodes to be
the weight of the minimum weight path between them.

HST metric spaces are particularly well-behaved when working with the
Wasserstein distance since the Wasserstein distance on a HST has a simple
closed form. A distribution $P$ on the the underlying metric space in a HST
induces a function $\mathfrak{G}_{P}$ on the nodes of the tree where the value
of a node $\nu$ is given by the weight in $P$ of the leaf nodes in the subtree
rooted at $\nu$. For every level $\ell\in[0pt]$ of the tree, let $P_{\ell}$ be
the distribution induced on the nodes at level $\ell$ where the probability of
node $\nu$ is $\mathfrak{G}_{P}(\nu)$. Thus $P_{\ell}$ is a discrete
distribution on a domain of size $N_{\ell}$, where $N_{\ell}$ is the number of
nodes in level $\ell$ of the tree.

######  Lemma 1.5 (Closed form Wasserstein distance formula).

Given two distributions $P$ and $Q$ defined in a HST metric space, the
Wasserstein distance between $P$ and $Q$ has the closed formula:

$\mathcal{W}(P,Q)=\frac{1}{2}\sum_{\nu}r_{\nu}|\mathfrak{G}_{P}(\nu)-\mathfrak{G}_{Q}(\nu)|=\sum_{\ell}r_{\ell}\text{\rm
TV}(P_{\ell},Q_{\ell}),$

where $r_{\nu}$ is the weight of the edge connecting $\nu$ to its parent, and
the sum is over all nodes in the tree.

We will call a node $\nu$ _$\alpha$ -active_ under the distribution $P$ if
$\mathfrak{G}_{P}(\nu)\geq\alpha$. Let $\gamma_{{P}}\left({\alpha}\right)$ be
the set of $\alpha$-active nodes under $P$ and
$\gamma_{{P_{\ell}}}\left({\alpha}\right)$ be the set of $\alpha$-active nodes
at level $\ell$. Then there exists an algorithm $\mathcal{A}$ such that given
a distribution $P$, $\varepsilon>0$, and $n\in\mathbb{N}$,

$\mathcal{R}_{\mathcal{A},n}(P)=\tilde{O}\left(\max_{\ell}r_{\ell}\sum_{x\in[N_{\ell}]}\min\left\\{P_{\ell}(x)(1-P_{\ell}(x)),\sqrt{\frac{P_{\ell}(x)(1-P_{\ell}(x))}{n}}\right\\}+\sum_{x\notin\gamma_{{P_{\ell}}}\left({2\kappa}\right)}P_{\ell}(x)+(|\gamma_{{P_{\ell}}}\left({2\kappa}\right)|-1)\kappa\right),$

where the max is over all the levels of the tree and
$\kappa=\Theta(\frac{\log(n)}{\varepsilon n})$. Further, this bound matches
(up to logarithmic factors) the lower bound $\min_{\varepsilon\text{-DP
}\mathcal{A}^{\prime}}\sup_{P^{\prime}:D_{\infty}(P,P^{\prime})\leq\ln
2}\operatorname*{\mathbb{E}}[\mathcal{W}(\mathcal{A}^{\prime}(\hat{P^{\prime}}_{n^{\prime}}),P^{\prime})]$
where $n^{\prime}\approx\frac{n}{\operatorname{polylog}n}$. The error rate
$\mathcal{R}_{\mathcal{A},n}$ does indeed adapt to easy instances as we
expected. The error decomposes into three components. The first component is
the non-private sampling error; the error that would occur even if privacy was
not required. The second component indicates that we can not privately
estimate the value of nodes that have probability less than $\approx
1/(\varepsilon n)$. The third component is the error due to privacy on the
active nodes. If $P$ is highly concentrated then we expect most nodes to
either be $\kappa$-active or have weight 0, so the first two terms in
$\mathcal{R}_{\mathcal{N},n,\varepsilon}(P)$ are small. There should also be
few active nodes, making the last term also small. Conversely, if $P$ has a
large region of low density then we expect a large number of inactive nodes,
as well as non-zero inactive nodes that are at higher levels of the tree and
hence contribute more to the final term. Thus, in distributions with high
dispersion we expect the right hand side to be large.

Upper Bounds: As in the one-dimensional setting, we want to restrict to only
privately estimating the density at a small number ($\approx\varepsilon n$) of
points. While we could try to mimic the one-dimensional solution by privately
estimating a solution to the $\varepsilon n$-median problem, it’s not clear
how to prove such an approach is instance-optimal. It turns out that a simpler
solution more amenable to analysis will suffice. Our algorithm has two stages;
first we attempt to find the set of $\kappa$-active nodes, then we estimate
the weight of these active nodes. Since these nodes have weight greater than
$\frac{\log(n)}{\varepsilon n}$, we can privately estimate them to within
constant multiplicative error. Any nodes that are not detected as active, are
initially ascribed a weight of 0. The error due to not estimating the non-
active nodes is absorbed into the third error term. The final step is to
project the noisy density function into the space of distributions on the
underlying metric space. The error of the upper bound algorithm is summed over
all levels of the tree, although since the depth of the tree is logarithmic in
the size of the metric space, this is within a logarithmic factor of the
maximum over the levels.

Lower Bound: We first observe that in order to estimate the distribution well
in Wasserstein distance, an algorithm must estimate each level of the tree
well in TV distance. This is derived from Lemma 1.5. This allows us to reduce
to the problem of lower bounding the error of density estimation of discrete
distributions in TV distance. The main tool we use is a differentially private
version of Assouad’s method. Similar to how the technique in the previous
section allowed us to relate lower bounding estimation rates to simple
hypothesis testing, Assouad’s lemma allows us to relate lower bounding
estimation rates to multiple hypothesis testing. Note that unlike the
technique in the previous section, Assouad’s lemma allows us to prove lower
bounds on the expected error, rather than lower bounds on high probability
error bounds. It involves constructing nets of distributions in
$\mathcal{N}(P)$ that are pairwise far in the relevant metric of interest
(which for us in the TV distance) but the multiple hypothesis testing problem
between the distributions is sufficiently hard. For proving the third term
belongs in the lower bound, the standard statement of DP Assouad’s lemma
[ASZ21] suffices, where one builds a set of distributions indexed by a
hypercube. For the first and second terms, we need to slightly generalise the
statement to allow for sets of distributions indexed by a product of
hypercubes. We use the approximate DP version of DP Assouad’s so while our
upper bounds are for pure differential privacy, our lower bounds hold for both
pure and approximate differential privacy.

Let us start with the third term. Suppose the number of active nodes is even
(a small tweak is made if there is an odd number of active nodes). We pair up
the active nodes and index each pair by a coordinate of the hypercube. For
each corner of the hypercube, $(u^{0},u^{1},\cdots,u^{k})\in\\{\pm 1\\}^{k}$,
for each coordinate $j\in[k]$, if $u^{j}=+1$, we move $\tilde{O}(\kappa)$ mass
from one node in the $j$th pair to the other node. If $u^{j}=-1$ then we leave
the $j$th pair of nodes alone. Since each active node has mass $>\kappa$, it’s
clear that each resulting distribution belongs in $\mathcal{N}(P)$. We can
also show that these distributions form a sufficiently hard multiple
hypothesis testing problem. By DP Assouad’s (Lemma 5.8), this allows us to
lower bound the estimation error by $\Omega(k\kappa)$, which is within
$\tilde{\Omega}$ of the third term when the number of active nodes is $\geq
2$. We treat the case where there is a single active node separately.

For the second term, we want to pair up the inactive nodes in a similar manner
and move half their mass from one node to the other. However, since we want to
remain within $\mathcal{N}(P)$, we can’t pair any two inactive nodes together.
Thus, we divide the inactive nodes into _scales_ , where nodes within a
certain scale all have weight within a multiplicative factor of two. We then
pair up nodes within each scale and have a different hypercube for each scale.
Again, it’s clear that these distributions are all in $\mathcal{N}(P)$ and we
can show that these distributions form a sufficiently hard multiple hypothesis
testing problem. The proof for the first term follows similarly.

##  2 Preliminaries

For all distributions $P$, we will use $f_{P}$ to denote the density of $P$
(when it exists) and $F_{P}$ to denote the cumulative distribution function of
$P$. Given a space $\mathcal{X}$, let $\Delta(\mathcal{X})$ be the set of
distributions on the space $\mathcal{X}$. Given a logical statement $a$, let
$\chi_{a}=0$ if $a$ is false and 1 if $a$ is true. For example, $\chi_{0=0}=1$
and $\chi_{0=1}=0$.

A number of distances between distributions are important in this work. We
start by defining the infinity divergence, which is important in the notion of
instance optimality we use.

######  Definition 2.1 ($D_{\infty}$-divergence).

Given two distributions $P$ and $Q$ with the same support, the $\infty$-Rényi
divergence
$D_{\infty}(P,Q)=\ln\sup_{t}\max\left\\{\frac{P(t)}{Q(t)},\frac{Q(t)}{P(t)}\right\\}$,
if $P$ and $Q$ are discrete, and
$D_{\infty}(P,Q)=\ln\sup_{t}\max\left\\{\frac{f_{P}(t)}{f_{Q}(t)},\frac{f_{Q}(t)}{f_{P}(t)}\right\\}$,
if $P$ and $Q$ are continuous distributions on $\mathbb{R}$, and have density
functions. If $P$ and $Q$ don’t have the same support, then
$D_{\infty}(P,Q)=\infty$.

We will use $\text{\rm KL}(P,Q)$ to denote the KL-divergence, $H^{2}(P,Q)$ to
denote the squared Hellinger divergence and $\text{\rm TV}(P,Q)$ to denote the
total variation distance, defined later. These metrics are defined in Appendix
A

##### Wasserstein Distance:

The error metric that we use to judge our performance on the density
estimation task is $1$-Wasserstein distance (that we will call just
Wasserstein distance where it is clear from context). In this subsection, we
define Wasserstein distance.

######  Definition 2.2.

For any separable metric space $(E,D)$, let $P,Q$ represent Borel measures on
$E$. Then, the $1$-Wasserstein distance between $P,Q$ is defined as

$\mathcal{W}(P,Q)=\inf_{\pi}\int_{E}\int_{E}D(t,t_{0})\pi(x,x_{0}),$

where the infimum is over all measures $\pi$ on the product space $E\times E$
with marginals $P$ and $Q$ respectively.

Finally, for one dimensional real spaces where the metric of interest is
$\ell_{1}$ norm, we will use the following equivalent formulation of
Wasserstein distance extensively.

######  Lemma 2.3 (Wasserstein formula over $\mathbb{R}$).

Let $P,Q$ represent probability distributions on $\mathbb{R}$ with finite
expectation. Then, the $1$-Wasserstein distance between $P,Q$ is equal to

$\mathcal{W}(P,Q)=\int_{\infty}^{\infty}|F_{P}(t)-F_{Q}(t)|dt,$

where the $F(\cdot)$ represents the cumulative distribution function.

Given an metric space $\mathcal{X}$, the Wasserstein metric is a well-defined
metric on the set of the probability distributions over $\mathcal{X}$.

###  2.1 Differential Privacy

We start by defining the Hamming distance between datasets.

######  Definition 2.4 (Hamming Distance).

For any $n>0$ and two datasets ${\bf x}=(x_{1},\dots,x_{n}),{\bf
x}^{\prime}=(x^{\prime}_{1},\dots,x^{\prime}_{n})\in\mathcal{X}^{n}$, the
Hamming distance $d_{Ham}({\bf x},{\bf x^{\prime}})$ between the two datasets
is defined as the number of entries that they disagree on, i.e.
$\sum_{i=1}^{n}\mathrm{1}[x_{i}\neq x^{\prime}_{i}]$.

Next, we define differential privacy.

######  Definition 2.5 (Differential Privacy [DMNS17, DKM+06]).

A randomized algorithm $\mathcal{A}:\mathcal{X}^{n}\rightarrow\mathcal{Y}$ is
said to be $(\varepsilon,\delta)$-differentially private if for every pair of
datasets ${\bf x},{\bf x}^{\prime}\in\mathcal{X}^{n}$ such that $d_{ham}({\bf
x},{\bf x}^{\prime})\leq 1$ and for all subsets $Y\subseteq\mathcal{Y}$,

$\Pr[\mathcal{A}({\bf x})\in Y]\leq e^{\varepsilon}\cdot\Pr[\mathcal{A}({\bf
x}^{\prime})\in Y]+\delta.$

Differential privacy satisfies two important properties that we will utilise:
it is closed under post-processing and composes naturally. More information
about these properties and an example of a basic differentially private
algorithm are given in Appendix A.

##  3 On Instance Optimality

In this section, we discuss the notion of instance optimality, and argue that
it provides a useful benchmark that captures the idea of _going beyond the
worst case_. The notion of instance optimality we propose can be see as a
generalisation of the hardest one-dimensional subproblem, or hardest local
alternative introduced by [CL15]. Suppose we have a family of distributions
$\mathcal{P}\subset\Delta(\mathcal{X})$ on a space $\mathcal{X}$ and our goal
is to learn the parameter $\theta:\mathcal{P}\to\mathcal{M}$ where
$\mathcal{M}$ is a metric space with metric $d$. Given an estimation algorithm
$\mathcal{A}:\mathcal{X}^{n}\to\mathcal{M}$, we can define the estimation
rate333We note that while the estimation rate here is defined in expectation,
we will sometimes show results (for e.g. in the one-dimensional case) where
estimation rate is defined with probability at least $0.75$ over the
randomness of the algorithm and the data; see Equation 1. of $\mathcal{A}$ to
be the function $\mathcal{R}_{\mathcal{A},n}:\mathcal{P}\to\mathbb{R}_{+}$
where

$\mathcal{R}_{\mathcal{A},n}(P)=\mathbb{E}_{D\sim
P^{n}}[d(\theta(P),\mathcal{A}(D))].$

Since the estimation rate is a function of the distribution $P$, the
estimation rate of an algorithm may be lower at “easy” distributions and
larger at “harder” distributions. As a classic example, consider the
estimation rate of Bernoulli parameter estimation where $\mathcal{A}$ simply
outputs the empirical mean. Then
$\mathcal{R}_{\mathcal{A},n}(\texttt{Ber}(p))=\min\\{p(1-p),\sqrt{p(1-p)/n}\\}$,
so this algorithm performs better when the Bernoulli parameter is close to 0
or 1, and has it’s worst case error when $p=1/2$.

Cai and Low [CL15] proposed three desiderata that a target estimation rate
$\mathcal{R}_{n}:\Delta(\mathcal{X})\to\mathbb{R}_{+}$ should satisfy in order
to be a meaningful benchmark;

  1. 1.

$\mathcal{R}_{n}(P)$ varies significantly across $\mathcal{P}$

  2. 2.

$\mathcal{R}_{n}$ is an achievable estimation rate; there exists an algorithm
$\mathcal{A}$ and constant $\alpha$ such that
$\mathcal{R}_{\mathcal{A},n}(P)\leq\alpha\mathcal{R}_{n}(P)$ for all
$P\in\mathcal{P}$

  3. 3.

Outperforming the benchmark $\mathcal{R}_{n}$ at one distribution leads to
worse performance at another distribution.

In this section we will discuss the definition of instance optimality we will
use in this work by defining the target estimation rate that will serve as our
benchmark estimation rate. The main theorems of this paper establish that our
chosen benchmark achieves desiderata 1 and 2 above. It is not immediately
obvious that desiderata 3 holds. We will show in Section 3.2, through the
introduction of a related notion of instance optimality which we call _local
minimality_ , that desiderata 3 holds in many important settings, including
the problem studied in this paper.

###  3.1 Local Estimation Rates

We will start by defining a target estimation rate. We’ll say an algorithm is
$\alpha$-instance optimal if it uniformly achieves this target estimation rate
up to a multiplicative $\alpha$ factor. For each distribution
$P\in\mathcal{P}$, we define a neighbourhood $\mathcal{N}(P)$.

######  Definition 3.1.

Given a function $\mathcal{N}:\mathcal{P}\to\mathfrak{P}(\mathcal{P})$, where
$\mathfrak{P}(\mathcal{P})$ is the power set of $\mathcal{P}$, we define the
optimal estimation rate with respect to $\mathcal{N}$ to be:

$\mathcal{R}_{\mathcal{N},n}(P)=\min_{\mathcal{A}}\sup_{Q\in\mathcal{N}(P)}\mathcal{R}_{\mathcal{A},n}(Q).$
(2)

An algorithm $\mathcal{A}$ is $\alpha$-instance optimal with respect to
$\mathcal{N}$ if for all $P\in\mathcal{P}$,

$\mathcal{R}_{\mathcal{A},n}(P)\leq\alpha\mathcal{R}_{\mathcal{N},n}(P)$

If an algorithm $\mathcal{A}$ uniformly achieves the optimal estimation rate
wrt a function $\mathcal{N}$, then this implies that for all distributions
$P$, the error of the algorithm $\mathcal{A}$ on $P$ is competitive with an
algorithm that is told the additional information that the distribution is in
$\mathcal{N}(P)$. Given a function $\mathcal{N}$, it is possible that there
does not exist an algorithm that uniformly achieves
$\mathcal{R}_{\mathcal{N},n}$. For example, as discussed in the introduction,
if $\mathcal{N}(P)=\\{P\\}$, then $\mathcal{R}_{\mathcal{N},n}$ is not
uniformly achievable. Conversely, if $\mathcal{N}(P)$ is not chosen carefully,
then the estimation rate $\mathcal{R}_{\mathcal{N},n}$ may not define a
meaningful benchmark; e.g. an estimation rate that adapts to easy instances.

A different formalization may be more probabilistic: the algorithm designer
may have in mind a distribution $\mathcal{D}$ over distributions that they
care about, and their objective may be to minimize
$\operatorname*{\mathbb{E}}_{P\sim\mathcal{D}}[\mathcal{R}_{\mathcal{A},n}(P)]$.
Suppose that for the $\mathcal{A}^{\star}$ chosen by the algorithm designer,
and for our neighborhood map $\mathcal{N}$, the function
$\mathcal{R}_{\mathcal{A}^{\star},n}(P)$ does not vary too much over
$\mathcal{N}(P)$ on average. Formally, let

$disc_{\mathcal{A}^{\star}}^{\mathcal{N}}(P)=\sup_{P^{\prime}\in\mathcal{N}(P)}(\mathcal{R}_{\mathcal{A}^{\star},n}(P^{\prime})-\mathcal{R}_{\mathcal{A}^{\star},n}(P)$

and let
$\overline{disc}_{\mathcal{A}^{\star}}^{\mathcal{N}}(P)=\operatorname*{\mathbb{E}}_{P\sim\mathcal{D}}[disc_{\mathcal{A}^{\star}}^{\mathcal{N}}(P)]$.
Then for any algorithm $\mathcal{A}$ that is $\alpha$-instance optimal with
respect to $\mathcal{N}$, we can write

$\displaystyle\operatorname*{\mathbb{E}}_{P\sim\mathcal{D}}[\mathcal{R}_{\mathcal{A},n}(P)]$
$\displaystyle\leq\alpha\cdot\operatorname*{\mathbb{E}}_{P\sim\mathcal{D}}[\sup_{P^{\prime}\in\mathcal{N}(P)}\mathcal{R}_{\mathcal{A}^{\star},n}(P^{\prime})]$
$\displaystyle=\alpha\cdot\operatorname*{\mathbb{E}}_{P\sim\mathcal{D}}[\mathcal{R}_{\mathcal{A}^{\star},n}(P)+disc_{\mathcal{A}^{\star}}^{\mathcal{N}}(P))]$
$\displaystyle=\alpha\cdot\operatorname*{\mathbb{E}}_{P\sim\mathcal{D}}[\mathcal{R}_{\mathcal{A}^{\star},n}(P)]+\alpha\cdot\operatorname*{\mathbb{E}}_{P\sim\mathcal{D}}[disc_{\mathcal{A}^{\star}}^{\mathcal{N}}(P))]$
$\displaystyle=\alpha\cdot\left(\operatorname*{\mathbb{E}}_{P\sim\mathcal{D}}[\mathcal{R}_{\mathcal{A}^{\star},n}(P)]+\overline{disc}_{\mathcal{A}^{\star}}^{\mathcal{N}}(P)\right).$

In other words, as long as the algorithm $\mathcal{A}^{\star}$’s performance
is relatively constant over $\mathcal{N}(P)$ on average over the distribution
of interest, the instance optimal algorithm (that is not tailored to
$\mathcal{D}$) is competitive with $\mathcal{A}^{\star}$. A similar result
holds for a multiplicative definition of $disc$.

This discussion can help guide the choice of the neighborhood function that is
appropriate for a particular application. In the case of density estimation in
the Wasserstein distance, we will define $\mathcal{N}(P)$ be a small
$D_{\infty}$ ball around $P$. We believe this captures the kind of domain
information an algorithm designer may have. E.g. one may have a small amount
of public data samples, in which case the posterior over distributions in a
$D_{\infty}$ ball will be relatively constant. If the algorithm designer’s
custom algorithm needs to do well for all distributions in this set, an
instance-optimal algorithm will be competitive with this custom algorithm.

Previous work in instance optimality has largely focused on two notions of
neighborhood. In [FLN01, ABC17, VV16, OS15], where the objects of interest are
discrete subsets with no a priori structure, it is natural to ask that the
algorithm work well for any permutation of the inputs. For example, if the
goal is to compute the set of maximal points from a 2-d point set, the
algorithm designer would typically want an algorithm that works well for any
permutation of the set of input points. In our setting where the points of
interest have a metric structure, this is not an appropriate notion. In fact,
even for the discrete case studied in Section 5.2.2, permutation invariance
cannot capture natural prior beliefs that may arise in practice. For example,
for power-law distributions that one often sees in private learning
applications [ZKM+20, CB22, CCD+23], a small number of samples are sufficient
to get a good estimate of the heavy bins, and rule out a large fraction of
permutations of the input space.

A second line of work arising from the statistics literature [CL15] has looked
at defining instance-optimality with respect to neighborhoods of size $2$.
While this approach has been very successful for many problems, we find it
inappropriate for density estimation (outside of density estimation on
$\mathbb{R})$ as neighborhoods of size two are too weak to capture the
difficulty of problems of interest. Even in the simple case of discrete
distributions, this neighborhood is provably insufficient to get instance-
optimality results with any $o(K)$ competitive ratio. Indeed, for any two
given distributions on $[K]$ with TV distance $\alpha$,
$\tilde{O}(\frac{1}{\alpha^{2}})$ samples suffice to distinguish them, whereas
learning a near uniform distribution on $K$ atoms requires $\Omega(K)$
samples. In the private setting, the need to use multiple distributions to
prove lower bounds is well-studied. Our approach shares this similarity of
using a multi-instance lower bounding argument with packing lower bounds in
privacy, and local Fano’s and Le Cam’s methods in statistics. Our work shows
that some of the same lower bounding techniques can be used to prove instance-
optimality results with respect to natural neighborhood maps, going well
beyond the the worst-case results those works prove.

In the special case of density estimation in the Wasserstein distance on
$\mathbb{R}$, instance optimality with respect to neighborhoods of size 2 is
achievable. In the standard version of this benchmark metric,
$\mathcal{N}(P)=\\{P,Q_{P}\\}$ where $Q_{P}$ can be _any_ distribution and is
chosen to maximise $\mathcal{R}_{\mathcal{N},n}(P)$. However, this notion may
not be an appropriate notion of instance optimality by itself. To see this,
consider a distribution $P$ supported on an interval $[a,b]$. Moving a small
amount of mass from one end of the interval to the other would create an
indistinguishable distribution that is far from $P$ in Wasserstein distance,
and a hypothesis testing argument can be used to show that the target
estimation rate defined above (for the hardest one-d sub problem) depends on
the interval size $b-a$. This implies that the adaptivity of algorithms to
support size of the distribution (crucial in Wasserstein estimation) is not
captured by this notion of instance optimality. Instead, we add a further
restriction to the definition to make it more appropriate for our setting; we
only consider distributions $Q$ that are in a small $D_{\infty}$ ball around
$P$ ($D_{\infty}(P,Q)\leq\ln 2)$, and ask that an algorithm is competitive
with an algorithm that is told the additional information that $P\in\\{P,Q\\}$
(in the worst case over distributions $Q$ that are in this $D_{\infty}$ ball).
That is, we define the benchmark estimation rate to be

$\mathcal{R}_{{\rm loc},n}(P)=\sup_{Q:D_{\infty}(P,Q)\leq\ln
2}\min_{\mathcal{A}}\\{\mathcal{R}_{\mathcal{A},n}(Q),\mathcal{R}_{\mathcal{A},n}(P)\\}.$
(3)

Note that all such distributions $Q$ have the same support as $P$, which
allows us to capture the adaptivity of algorithms to the support size of the
distribution. Specifically, we define the following target estimation rate in
the one-dimensional setting. In the case of estimating distributions on a
bounded subset of $\mathbb{R}$, we will show that this error rate is
achievable, up to logarithmic factors.

We also note that our notion of instance optimality more naturally captures
the accuracy of algorithms even for basic tasks. Note that for the Bernoulli
case, our technique achieves a bound of
$\frac{\sqrt{p(1-p)}}{\sqrt{n}}+\min\\{p,1-p,\frac{1}{\varepsilon n}\\}$ which
also appear to be better than the instance-optimal lower bounds in [MSU22],
which take the form $\frac{\sqrt{p(1-p)}}{\sqrt{n}}+\frac{1}{\varepsilon n}$.
This apparent contradiction can be explained by the the use in [MSU22] of the
hardest-one dimensional sub-problem to define the instance-optimal rate, i.e.,
$\mathcal{N}(P)$ is $\\{P,Q\\}$ for a worst-case Bernoulli $Q$. On the other
hand, the notion of instance-optimality we use would only consider Bernoullis
$Q$ such that $D_{\infty}(P,Q)\leq\ln 2$. When $p$ is close to $0$, the lower
bound in [MSU22] would on this instance consider $Q$ to be
$Bern(p+\frac{1}{\varepsilon n})$, which can have a large
$D_{\infty}$-distance from $P$, and so isn’t in the neighborhood used in our
notion of instance-optimality. Hence, the target rate one would obtain from
our definition is smaller when $p$ is close to $0$. Our algorithm can achieve
this improved rate, as it is likely to output $0$ as an estimate of $p$ in
this case, pushing small counts down to zero.

Recent differentially private algorithms such as those in [HLY21, DKSS23] have
shown instance-optimality for problems such as mean estimation. Relatedly,
other works have designed algorithms that adapt to the local/smooth/deletion
sensitivity of the underlying function. An instance in these works in a
dataset rather than a distribution, and it is not clear how to extend the
corresponding notion of neighborhood to our setting. Our neighborhood notion
perhaps comes closest to the deletion neighborhoods considered in some of
these works.

Finally, we remark that while we have stated our results as being competitive
with the worst-case instance in $\mathcal{N}(P)$, they apply for the average
case over a specific distribution over $\mathcal{N}(P)$. Since that specific
distribution is adversarial, we don’t view this version as more natural than
the worst case.

Given that we are focusing on private estimation, we will use use
$\mathcal{R}_{{\rm loc},n,\varepsilon}$ to denote the version of Eqn 3 where
the minimum is taken over all $\varepsilon$-DP mechanisms, and
$\mathcal{R}_{\mathcal{N},n,\varepsilon}$ to define the optimal
$\varepsilon$-DP estimation rate, i.e. Eqn 2 where the minimum is taken over
all $\varepsilon$-DP mechanisms.

###  3.2 Locally Minimal Algorithms

In this section we address the third desiderata of [CL15]. An important
concept in statistics is that of efficiency of an estimator, which informally
compares the rate of convergence of the estimator with a benchmark that in
general is not beatable. This idea has been used to argue that for some
fundamental estimation problems, the Maximum Likelihood Estimator (MLE) is the
best possible. Hodge showed an example of a superefficient estimator that is
asymptotically as good as the MLE everywhere, but beats the MLE on a certain
set of inputs. The statistics community has argued in multiple ways that these
superefficient estimators do not limit our ability to argue that MLE is
“optimal”. We refer the reader to [vdV97, Wol65, Vov09] for a discussion of
superefficiency. One of the more compelling arguments here is a result saying
that the set of points where superefficiency is achieved has Lebesgue measure
zero. This in particular implies that in a small neighborhood around any
point, there is a point (in fact many points) where the superefficient
estimator does no better than the MLE. In the partial order on estimators, the
MLE is thus minimal and this is true even when looking at the performance of
the estimator only on a small neighborhood around a given point.

This motivates a slightly different notion capturing the goodness of the
algorithm locally.

######  Definition 3.2.

Let $\mathcal{M}$ be a class of algorithms. We say that an algorithm
$\mathcal{A}$ is $\alpha$-locally minimal with respect to a neighborhood map
$\mathcal{N}$, if for all instance $P$, and all
$\mathcal{A}^{\prime}\in\mathcal{M}$, there is a $Q\in\mathcal{N}(P)$ such
that
$\mathcal{R}_{\mathcal{A},n}(Q)\leq\alpha\cdot\mathcal{R}_{\mathcal{A}^{\prime},n}(Q)$.

In words, local minimality says that for any other $\mathcal{A}^{\prime}$, the
algorithm $\mathcal{A}$ is competitive with $\mathcal{A}^{\prime}$ for some
instance in the neighborhood of $P$. Put differently, no
$\mathcal{A}^{\prime}$ can be uniformly much better than $\mathcal{A}$ on the
neighborhood, even one that knows $P$.

We show that in general, this notion is incomparable to our notion of instance
optimality. Nevertheless, under reasonable assumptions, the two notions are
closely related.

######  Example 3.3 (Local Minimality $\not\Rightarrow$ Instance Optimality).

Consider a pair of instances $\\{P,Q\\}$ with
$\mathcal{N}(P)=\mathcal{N}(Q)=\\{P,Q\\}$. Let $\mathcal{M}$ contain two
algorithms $\mathcal{A}$, and $\mathcal{A}^{\star}$ with

$\displaystyle\mathcal{R}_{\mathcal{A},n}(P)=1;\;\;\;\;\;\;$
$\displaystyle\mathcal{R}_{\mathcal{A},n}(Q)=0;$
$\displaystyle\mathcal{R}_{\mathcal{A}^{\star},n}(P)=0;\;\;\;\;\;\;$
$\displaystyle\mathcal{R}_{\mathcal{A}^{\star},n}(Q)=0;$

Then one can verify that $\mathcal{A}$ is ($1$-)locally minimal in
$\mathcal{M}$. However, it is not $\alpha$-instance optimal for any finite
$\alpha$ as it fails to satisfy the definition at $P$.

######  Example 3.4 (Instance Optimality $\not\Rightarrow$ Local Minimality).

Consider a set of instances $\\{P_{1},P_{2},P_{3}\\}$ with
$\mathcal{N}(P_{1})=\\{P_{1},P_{2}\\},\mathcal{N}(P_{2})=\\{P_{1},P_{2},P_{3}\\},\mathcal{N}(P_{3})=\\{P_{2},P_{3}\\}$.
Let $\mathcal{M}$ contain algorithms $\mathcal{A},\mathcal{A}^{\star}$ with

$\displaystyle\mathcal{R}_{\mathcal{A}^{\star},n}(P_{1})=1;\;\;\;\;\;\;$
$\displaystyle\mathcal{R}_{\mathcal{A},n}(P_{1})=2\alpha;$
$\displaystyle\mathcal{R}_{\mathcal{A}^{\star},n}(P_{2})=2\alpha;\;\;\;\;\;\;$
$\displaystyle\mathcal{R}_{\mathcal{A},n}(P_{2})=4\alpha^{2};$
$\displaystyle\mathcal{R}_{\mathcal{A}^{\star},n}(P_{3})=4\alpha^{2};\;\;\;\;\;\;$
$\displaystyle\mathcal{R}_{\mathcal{A},n}(P_{3})=4\alpha^{2}.$

Then one can verify that $\mathcal{A}$ is ($1$-)instance optimal in
$\mathcal{M}$. However, it is not $\alpha$-locally minimal at $P_{1}$.

Under smoothness assumptions on $\mathcal{A}$ with respect to $\mathcal{N}$,
one can argue that the two notions are essentially equivalent.

######  Proposition 3.5.

Let $\mathcal{A}$ be such that for all instances $P$ and for all
$Q\in\mathcal{N}(P)$,
$\mathcal{R}_{\mathcal{A},n}(Q)\leq\beta\cdot\mathcal{R}_{\mathcal{A},n}(P)$.
Further, suppose that $\mathcal{N}(P)$ is compact for any $P$. If
$\mathcal{A}$ is $\alpha$-instance optimal in $\mathcal{M}$ with respect to
$\mathcal{N}$, then it is $\alpha\beta$-locally minimal.

###### Proof.

Let $P$ be an instance and let $\mathcal{A}^{\prime}$ be a competing
algorithm. By definition of $\alpha$-instance optimality,

$\displaystyle\mathcal{R}_{\mathcal{A},n}(P)$
$\displaystyle\leq\alpha\cdot\sup_{Q\in\mathcal{N}(P)}\mathcal{R}_{\mathcal{A}^{\prime},n}(Q).$

By compactness, this implies that there is a $Q$ achieving the supremum. In
other words, there exists $Q^{\star}\in\mathcal{N}(P)$ such that

$\displaystyle\mathcal{R}_{\mathcal{A},n}(P)$
$\displaystyle\leq\alpha\cdot\mathcal{R}_{\mathcal{A}^{\prime},n}(Q^{\star}).$

Since $Q^{\star}\in\mathcal{N}(P)$, our smoothness assumption implies that

$\displaystyle\mathcal{R}_{\mathcal{A},n}(Q^{\star})$
$\displaystyle\leq\beta\cdot\mathcal{R}_{\mathcal{A},n}(P).$

Combining the last two inequalities, this $Q^{\star}$ satisfies

$\displaystyle\mathcal{R}_{\mathcal{A},n}(Q^{\star})$
$\displaystyle\leq\alpha\beta\cdot\mathcal{R}_{\mathcal{A}^{\prime},n}(Q^{\star}).$

Since $P$ and $\mathcal{A}^{\prime}$ were arbitrary, this implies that
$\mathcal{A}$ is $\alpha\beta$-locally minimal. ∎

######  Proposition 3.6.

Let $\mathcal{A}$ be such that for all instances $P$ and for all
$Q\in\mathcal{N}(P)$,
$\mathcal{R}_{\mathcal{A},n}(Q)\geq\beta^{-1}\cdot\mathcal{R}_{\mathcal{A},n}(P)$.
If $\mathcal{A}$ is $\alpha$-locally minimal in $\mathcal{M}$ with respect to
$\mathcal{N}$, then it is $\alpha\beta$-instance optimal.

###### Proof.

Let $P$ be an instance and let $\mathcal{A}^{\prime}$ be a competing
algorithm. By definition of $\alpha$-local minimality, there is a
$Q^{\star}\in\mathcal{N}(P)$ such that

$\displaystyle\mathcal{R}_{\mathcal{A},n}(Q^{\star})$
$\displaystyle\leq\alpha\cdot\mathcal{R}_{\mathcal{A}^{\prime},n}(Q^{\star}).$

Since $Q^{\star}\in\mathcal{N}(P)$, our smoothness assumption implies that

$\displaystyle\mathcal{R}_{\mathcal{A},n}(P)$
$\displaystyle\leq\beta\cdot\mathcal{R}_{\mathcal{A},n}(Q^{\star}).$

Combining the last two inequalities, this $Q^{\star}$ satisfies

$\displaystyle\mathcal{R}_{\mathcal{A},n}(P)$
$\displaystyle\leq\alpha\beta\cdot\mathcal{R}_{\mathcal{A}^{\prime},n}(Q^{\star})$
$\displaystyle\leq\alpha\beta\cdot\sup_{Q\in\mathcal{N}(P)}cost(\mathcal{A}^{\prime}(Q),Q).$

Since $P$ and $\mathcal{A}^{\prime}$ were arbitrary, this implies that
$\mathcal{A}$ is $\alpha\beta$-instance optimal. ∎

A similar pair of results hold when the comparator algorithm
$\mathcal{A}^{\prime}$ is smooth with respect to the neighborhood map.

###  3.3 Relaxed Definitions

We finish by noting relaxations of the above definitions that share the same
semantic meaning (our algorithms will achieve these relaxed notions).

######  Definition 3.7.

Given a function $\mathcal{N}:\mathcal{P}\to\mathfrak{P}(\mathcal{P})$, where
$\mathfrak{P}(\mathcal{P})$ is the power set of $\mathcal{P}$, we define the
optimal estimation rate with respect to $\mathcal{N}$ to be:

$\mathcal{R}_{\mathcal{N},n}(P)=\min_{\mathcal{A}}\sup_{Q\in\mathcal{N}(P)}\mathcal{R}_{\mathcal{A},n}(Q).$
(4)

An algorithm $\mathcal{A}$ is $(\alpha,\beta,\gamma)$-instance optimal with
respect to $\mathcal{N}$ if for all $P\in\mathcal{P}$,

$\mathcal{R}_{\mathcal{A},n}(P)\leq\alpha\mathcal{R}_{\mathcal{N},\beta
n}(P)+\gamma$

######  Definition 3.8.

Let $\mathcal{M}$ be a class of algorithms. We say that an algorithm
$\mathcal{A}$ is $(\alpha,\beta,\gamma)$-locally minimal with respect to a
neighborhood map $\mathcal{N}$, if for all instance $P$, and all
$\mathcal{A}^{\prime}\in\mathcal{M}$, there is a $Q\in\mathcal{N}(P)$ such
that
$\mathcal{R}_{\mathcal{A},n}(Q)\leq\alpha\cdot\mathcal{R}_{\mathcal{A}^{\prime},\beta
n}(Q)+\gamma$.

Note that we think of $\beta\in(0,1]$ and $\gamma$ as non-negative. The reason
these are relaxed definitions is because we allow for an additive
approximation factor in addition to a multiplicative factor, and also compare
to a benchmark rate that depends on a potentially smaller number of samples
(and is hence easier to achieve). The original definition of instance
optimality (Definition 3.1) can be obtained by setting $\beta=1$ and
$\gamma=0$.

In our work, for most settings of interest, we roughly achieve $\beta=1/(\log
n)^{O(1)}$ and $\gamma$ to be an arbitrarily small polynomial in the inverse
of the number of samples $1/n$ at a $\log(1/\gamma)$ cost to the
multiplicative factor. We don’t view this as a significant issue since we
expect the benchmark rate with $\tilde{O}(n/\log n)$ samples to behave
asymptotically similarly to that with $n$ samples in most cases. We leave it
as an open question as to whether the original definition of instance
optimality can be achieved.

##  4 Additional Related Work

##### Instance Optimality for Differentially Private Statistics:

Several recent works have focused on formulating and giving ‘instance optimal’
differentially private algorithms for various statistical tasks. The work of
McMillan, Smith and Ullman [MSU22] is most directly related to our work; they
gave locally minimax optimal algorithms for parameter estimation for one-
dimensional exponential families in the central model of differential privacy.
The work of Duchi and Ruan [DR18] also gives locally-minimax optimal
algorithms for various one-dimensional parameter estimation problems under the
stronger constraint of local differential privacy. The notion of local minimax
optimality both these papers use is based on the hardest one-dimensional sub-
problem described in Section 3.1. While our results for density estimation in
$\mathbb{R}^{1}$ satisfy this notion, they also satisfy a stronger notion
described in Section 3.1. Additionally, as discussed in [MSU22], this
definition is provably unsuitable for higher dimensions; we instead suggest a
looser definition of instance optimality that is more promising in higher
dimensions. More importantly, our paper is primarily focused on the non-
parametric setting, and hence our techniques are different than the ones used
in those papers, which focused primarily on parameter estimation.

##### Other Beyond Worse-Case Results in Central Differential Privacy:

Several additional works in the differential privacy literature study
algorithms with accuracy that varies with the input dataset. Nearly all of
them look at the empirical setting where we are concerned with the specific
input dataset, rather than a distribution it may be drawn from. While initial
algorithms in differential privacy added noise based on a worse case notion of
global sensitivity, these works give various algorithmic frameworks that help
develop algorithms with guarantees that adapt to the hardness of the input
dataset. These include algorithms based on smooth sensitivity [NRS07, BS19],
the propose-test-release framework [DL09, BA20], Lipschitz extensions [BBDS13,
KNRS13, CZ13, RS16], and sensitivity pre-processing [CD20]. However, none of
these works study a formal notion of instance optimality.

In contrast, some more recent work do study definitions of instance optimality
in the empirical setting. A work of Asi and Duchi [AD20] studies two notions
of instance optimality: one by comparing the performance of an algorithm on a
dataset against the performance of the best unbiased algorithm on that
dataset, and another based on an analogue of the ‘hardest one-dimensional sub-
problem’ for the empirical setting (they compare the performance of an
algorithm on a dataset with all benchmark algorithms that know that the input
dataset is either of two possible datasets but whose performance is evaluated
as the worse over the two datasets). They give a general mechanism known as
the inverse sensitivity mechanism that they show is nearly instance optimal
under these definitions for various problems such as median and mean
estimation. Our work is focused on population quantities as opposed to
empirical quantities—while these are related, they can be very different. For
example, as pointed out in McMillan, Smith and Ullman [MSU22], using the
inverse sensitivity mechanism in [AD20] to estimate the mean of a Gaussian (by
using a locally minimax optimal algorithm for empirical mean) will result in
infinite mean squared error, whereas other approaches that reason directly
about the population quantities can get much better error.

In [DKSS23] and [HLY21], different notions of instance optimality are defined.
Roughly, they compare the performance of an algorithm on a dataset with a
benchmark algorithm that knows the input dataset but whose performance is
evaluated as the worst-case performance over large subsets of the input
dataset. While the details of the definitions in these papers vary slightly,
both papers give instance-optimal algorithms for mean estimation under their
respective definitions. For one-dimensional distributions, our algorithmic
technique at a high level shares ideas with these algorithms—the algorithms in
their papers try to adapt to the range of values in the dataset, whereas we
try to adapt to the level of concentration of the distribution. However, the
details of how this is done and the associated analyses vary. Our algorithm
for general metric spaces uses different techniques. Our work differs from
these works in a few other prominent ways: firstly, they are primarily
concerned with estimating functionals of the underlying dataset, whereas we
are concerned with density estimation in Wasserstein distance—these are
problems with different output types and different error metrics. Finally, it
is not clear if notions such as subset-based instance optimality that are well
defined in the empirical setting transfer meaningfully to the distributional
setting.

##### Instance-Optimal Statistical Estimation without Privacy Constraints:

Donoho and Liu [DL91] formulated the notion of the ‘hardest one-dimensional
sub-problem’ as a way of capturing instance optimality for statistical
estimation and gave non-private instance optimal algorithms for some one-
dimensional parameter estimation problems. Cai and Low [CL15] formulated an
instance-optimality type definition for non-parameteric estimation problems.
Our results for Wasserstein density estimation over $\mathbb{R}$ use a
stronger version of this notion of instance optimality. In higher dimensions,
this notion is provably unachievable, and so we define a different notion.

The other line of work most related to ours is on instance-optimal learning of
discrete distributions [OS15, VV16, HO19]. In their setting, instance
optimality is defined by comparing the performance of an algorithm on a
discrete distribution $P$ to the minimax error of any algorithm on the class
of discrete distributions with probability vectors that are permutations of
the probability vector of $P$. We note that this notion is not well suited to
many metric spaces, because permutations may not preserve properties such as
concentration of the distribution, and hence this notion of instance
optimality may provide an overly pessimistic view of the performance of an
algorithm. Our notion of instance optimality (in terms of $D_{\infty}$
neighborhood) compares against algorithms with a different type of prior
knowledge- i.e., the location of where the distribution concentrates, and
approximate values of the probabilities at each point. We note that these are
technically incomparable, and may be useful in different settings. For
estimation in Wasserstein distance, knowledge of where the distribution is
concentrated could be very useful in algorithm design, and so comparing to
algorithms with this type of knowledge is more appropriate. See Section 3 for
more discussion.

Finally, there is another line of work on getting similar instance optimal
guarantees for other statistical problems [ADJ+11, ADJ+12, AJOS13b, AJOS13a].
For the closeness testing problem (given two sequences, determine if they are
produced by the same distribution, or different distributions), Acharya, Das,
Jafarpour, Orlitsky, Pan and Suresh [ADJ+11, ADJ+12] developed a test (without
any knowledge about the generating distributions) that achieves the same error
with $O(n^{3/2})$ samples that an optimal label-invariant test that knows the
distributions $p$ and $q$ would achieve with $n$ samples.

##### Other work on Differentially Private Statistics:

There is a lot of other work on private statistical estimation, and we survey
the most relevant parts of the literature here. There is a long line of work
on minimax parameter/distribution estimation on various parametric
distribution families: product distributions [BUV18, KLSU19, ASZ20, CWZ19,
Sin23], Gaussian, sub-Gaussian distributions (and more generally exponential
families) [KV18, KLSU19, AAK21, BGS+21, KMS+22b, KMS22a, HKM22, KMV22, AL22,
LKO22, TCK+22, HKMN23, AKT+23, BHS23, KDH23], mixtures of Gaussian
distributions [KSSU19, AAL23b, AAL23a], heavy-tailed distributions [KSU20,
Nar23], discrete distributions with finite support [DHS15, ASZ20],
distributions with finite covers [BKSW21] and more. This line of work focuses
on minimax guarantees in the parametric setting, i.e. optimizing the worst-
case error of an algorithm over the entire class of distributions. Our work,
on the other hand works in the non-parametric setting where we do not make
assumptions about the distribution the dataset is drawn from, but instead give
‘instance-optimal’ algorithms that adapt to the hardness of the distribution
the input dataset is drawn from.

There is also a line of work on differentially private CDF estimation [DNPR10,
CSS11, BNS16, BNSV15, ALMM19, KLM+20, CLN+23], and quantile estimation [KSS22,
GJK21, ASSU24]. Our algorithm for density estimation over $\mathbb{R}$ uses a
quantile estimation algorithm (based on a CDF estimator) as a subroutine.
Finally, there is a line of work on differentially private testing [ASZ17,
CDK17, CKM+19], and the work characterizing the sample complexity of simple
hypothesis tests forms an important part of our analysis of the instance-
optimal rate for distributions over $\mathbb{R}$.

##### Work on Estimation in Wasserstein Distance:

In addition to the recent works [BSV22, HVZ23] on private Wasserstein learning
on $[0,1]^{d}$, there is a plethora of works studying it in the non-private
setting.

One line of work studies the convergence in Wasserstein distance of the
empirical measure (on $n$ samples) to the true measure, as a function of the
measure and the number of samples $n$ [Dud69, DY95, CR12, DSS11, BG14, FG15,
BL19, WB19, Lei20, Fou23]. Some of the later works above can be viewed as
studying this problem from a beyond worst-case analysis viewpoint. They give
upper and lower bounds for the expected value of this quantity, in terms of
various notions of ‘dimension’ of the underlying measure, such as the covering
number of the support of the distribution, the upper and lower ‘Wasserstein
dimensions’ of the measure, and others. Our work shows that the empirical
measure, appropriately massaged, is approximately instance-optimal for density
estimation without privacy constraints (for the notions of instance optimality
we consider), and hence these works give us a handle on the instance-optimal
rate as a function of the distribution and sample size $n$. Some more recent
work studies minimax estimation in Wasserstein distance [SP19, NWB19], and
show that without additional assumptions on the distribution, the empirical
measure is minimax optimal. Our work extends this result to show that in the
general non-parametric setting, the empirical measure is also approximately
instance-optimal; to the best of our knowledge, instance optimal estimation in
Wasserstein distance (even without privacy constraints) has not been
previously studied.

##  5 Distribution Estimation on Hierarchically Separated Trees

Let us now turn to distribution estimation on arbitrary finite metric spaces.
We will use the fact that any metric on a finite space can be embedding in a
hierarchically separated tree (HST) metric to reduce the problem of density
estimation in Wasserstein distance on an arbitrary metric space to density
estimation in Wasserstein distance on an HST. In Section 5.2 we’ll
characterise the target estimation rate $\mathcal{R}_{\mathcal{N},n}$. In
Section 5.3, we’ll then provide an $\varepsilon$-DP algorithm and prove that
it achieves this target estimation rate up to logarithmic factors.

###  5.1 Preliminaries on Hierarchically Separated Trees

A key component of our proof strategy is the reduction to Hierarchically
Separated Trees (HSTs). HSTs are special class of tree metrics that are able
to embed arbitrary metric spaces with low distortion. They are particularly
well-behaved when working with the Wasserstein distance since the Wasserstein
distance on an HST has a simple closed form.

######  Definition 5.1 (Hierarchically Separated Tree).

A hierarchically separated tree (HST) is a rooted weighted tree such that the
edges between level $\ell$ and $\ell-1$ all have the same weight (denoted
$r_{\ell}$) and the weights are geometrically decreasing so
$r_{\ell+1}=(1/2)r_{\ell}$. Let $0pt$ be the depth of the tree.

An HST defines a metric on its leaf nodes by defining the distance between any
two leaf nodes to be the weight of the minimum weight path between the two
nodes. We will rely on two main facts about HSTs in this work.

######  Lemma 5.2 (Low distortion metric embeddings [FRT03]).

Let $(V,d)$ be a metric space with $M$ points. There exists a randomized,
polynomial time algorithm that produces an HST where the leaf nodes of the
tree correspond to the elements of the metric space and the induced tree
metric $d_{T}$ is such that for all $u,v\in V$

  * •

$d(u,v)\leq d_{T}(u,v)$

  * •

$\mathbb{E}[d_{T}(u,v)]\leq O(\log M)\cdot d(u,v)$

The depth of the HST is logarithmic in the size of the metric space, $0pt=\log
M$.

An immediate consequence of the $O(\log M)$ metric distortion in Lemma 5.2 is
that the Wasserstein distance in the original metric space is also preserved
up to a $O(\log M)$ factor in expectation. Thus, Lemma 5.2 allows us to
translate the problem of learning densities on an arbitrary metric space in
Wasserstein distance to learning densities in Wasserstein distance on an HST.
This is a useful tool since HST metrics are generally easier to work with and,
as we’ll see below, the Wasserstein distance is particularly well-behaved on
an HST. In order to use Lemma 5.2 to translate the problem of density
estimation on a bounded ball in $\mathbb{R}^{d}$ into density estimation on an
HST, one discretizes the metric, paying a small additive term.

######  Corollary 5.3.

Given $\alpha>0$, there is a probabilistic embedding $f$ of $[0,1]^{d}$ into
an HST such that for all $x,y\in[0,1]^{d}$:

  * •

$d(x,y)-\alpha\leq d_{T}(f(x),f(y))$

  * •

$\mathbb{E}[d_{T}(f(x),f(y)]\leq
O(d\cdot\log\frac{1}{\alpha}])\cdot(d(x,y)+\alpha)$

The distortion is logarithmic in $\frac{1}{\alpha}$, so taking $\alpha$ to be
polynomially small, one gets the distortion to be $O(d\log n)$. It is easy to
see that this implies that the Wasserstein distance is preserved in both
directions up to $O(d\log\frac{1}{\alpha}$, up to an $\alpha$ additive error.

A distribution $P$ on the the underlying metric space in an HST induces a
function $\mathfrak{G}_{P}$ on the nodes of the tree where the value of a node
$\nu$ is given by the weight in $P$ of the leaf nodes in the subtree rooted at
$\nu$. For every level $\ell\in[0pt]$ of the tree, let $P_{\ell}$ be the
distribution induced on the nodes at level $\ell$ where the probability of
node $\nu$ is $\mathfrak{G}_{P}(\nu)$. Thus $P_{\ell}$ is a discrete
distribution on a domain of size $N_{\ell}$, where $N_{\ell}$ is the number of
nodes in level $\ell$ of the tree.

######  Lemma 5.4 (Closed form Wasserstein distance formula).

Given two distributions $P$ and $Q$ defined on an HST metric space, the
Wasserstein distance between $P$ and $Q$ has the closed formula:

$\mathcal{W}(P,Q)=\frac{1}{2}\sum_{\nu}r_{\nu}|\mathfrak{G}_{P}(\nu)-\mathfrak{G}_{Q}(\nu)|=\sum_{\ell}r_{\ell}\text{\rm
TV}(P_{\ell},Q_{\ell}),$

where $r_{\nu}$ is the length of the edge connecting $\nu$ to its parent, and
the sum is over all nodes in the tree.

###  5.2 The Target Estimation Rate

Recall the definition of our neighbourhood.

$\mathcal{N}(P)=\\{Q\in\mathcal{P}\;|\;D_{\infty}(P,Q)\leq\ln 2\\}$

We will call a node $\nu$, _$\alpha$ -active node under the distribution $P$_
if the weight in $P$ of the sub-tree rooted at $\nu$ is greater than $\alpha$.
Let $\gamma_{{P}}\left({\alpha}\right)$ be the set of $\alpha$-active nodes
under $P$ and $\gamma_{{P_{\ell}}}\left({\alpha}\right)$ be the
$\alpha$-active nodes at level $\ell$.

######  Theorem 5.5.

Given a distribution $P$ on $[N]$, $\varepsilon>0$, $\delta\in[0,1]$, and
$n\in\mathbb{N}$, let $\kappa=\frac{1}{10\varepsilon
n}\min\\{W\left(\frac{0.45\varepsilon}{\delta}\right),0.6\\}$ where $W(x)$ is
the Lambert W function so $W(x)e^{W(x)}=x$, then

$\mathcal{R}_{\mathcal{N},n,\varepsilon}(P)=\Omega\left(\max_{\ell}r_{\ell}\sum_{x\in[N_{\ell}]}\min\left\\{P_{\ell}(x)(1-P_{\ell}(x)),\sqrt{\frac{P_{\ell}(x)(1-P_{\ell}(x))}{n}}\right\\}+\sum_{x\notin\gamma_{{P_{\ell}}}\left({2\kappa}\right)}P_{\ell}(x)+(|\gamma_{{P_{\ell}}}\left({2\kappa}\right)|-1)\kappa\right),$

where the max is over all the levels of the tree.

Note that $\kappa\approx\frac{1}{\varepsilon n}\min\\{\log(1/\delta),1\\}$ so
the dependence on $\varepsilon$ and $n$ in Theorem 5.5 matches the upper bound
in Theorem 5.13. The error rate $\mathcal{R}_{\mathcal{N},n,\varepsilon}$ does
indeed adapt to easy instances as we expected. The error decomposes into three
components. The first component is the non-private sampling error; the error
that would occur even if privacy was not required. The second component
indicates that we can not estimate the value of nodes that have probability
less than $1/(\varepsilon n)$. The third component is the error due to privacy
on the active nodes. If $P$ is highly concentrated then we expect most nodes
to either be $\frac{1}{\varepsilon n}$-active or have weight 0, so the first
two terms in $\mathcal{R}_{\mathcal{N},n,\varepsilon}(P)$ are small. There
should also be few active nodes, making the last term smaller as well.
Conversely, if $P$ has a large region of low density then we expect a large
number of inactive nodes, as well as non-zero inactive nodes that are at
higher levels of the tree and hence contribute more to the final term. Thus,
in distributions with high dispersion we expect the right hand side to be
large.

The proof of Theorem 4.1 will involve two main steps. First, we will reduce
the lower bound on the HST to a lower bound on a star metric, or equivalently
estimation of a discrete distribution in TV distance. We’ll then use a variant
of Assouad’s inequality to prove the lower bounds on estimating discrete
distributions in TV distance.

####  5.2.1 Reduction to Estimation in TV distance of Discrete Distributions

The key observation is that in order to estimate the distribution well in
Wasserstein distance, an algorithm must estimate each level of the tree well
in TV distance. Any estimate of $P$ also induces an estimate of $P_{\ell}$;
let $\hat{P}$ be an estimate of the distribution $P$ and $\hat{P}_{\ell}$ be
the induced estimate of the distribution at level $\ell$. Then for any
distribution $P$

$\mathcal{W}(P,\hat{P})=\sum_{\ell\in[0pt]}r_{\ell}TV(P_{\ell},\hat{P}_{\ell}).$

The following observation ensures that our notions of instance optimality in
both the Wasserstein metric and the per-level TV distance are compatible at
every level $\ell$.

######  Theorem 5.6.

For every level $\ell\in[0pt]$, define the neighborhood of $P_{\ell}$ as
$\mathcal{N}_{\ell}:\Delta([N_{\ell}])\to\mathfrak{P}(\Delta([N_{\ell}]))$ by
$\mathcal{N}_{\ell}(P_{\ell})=\\{Q_{\ell}\;|\;D_{\infty}(P_{\ell},Q_{\ell})\leq\ln
2\\}$. Then,

$\mathcal{R}_{\mathcal{N},n,\varepsilon}(P)\geq\max_{\ell\in[0pt]}r_{\ell}\cdot\mathcal{R}_{\mathcal{N}_{\ell},n,\varepsilon}(P_{\ell}),$

where the error of $P$ is measured in the Wasserstein distance and $P_{\ell}$
is measured in the TV distance.

Recall that $\mathcal{R}_{\mathcal{N}_{\ell},n,\varepsilon}(P_{\ell})$ is the
optimal estimation rate with respect to $\mathcal{N}_{\ell}$ where the error
is measured with respect to the total variation error. The proof of Theorem
5.6 can be found in Appendix C.

####  5.2.2 Characterizing Target Estimation Rate for Discrete Distributions

In light of Theorem 5.6, we will focus on characterizing the difficulty of
estimating the distribution at a single level of the tree for the remainder of
this section. Since this is fundamentally a statement about estimating
discrete distributions in TV distance, we will state everything in this
section in terms of general discrete distributions. Let $N\in\mathbb{N}$, and
let $P$ be a distribution on $[N]$. Define
$\mathcal{N}(P)=\\{Q\;|\;D_{\infty}(P,Q)\leq\ln 2\\}$. Our goal is to give a
lower bound for $\mathcal{R}_{f,n,\varepsilon}(P)$, where the metric is the TV
distance.

######  Theorem 5.7.

Given $\varepsilon>0$ and $\delta\in[0,1]$, let $\kappa=\frac{1}{10\varepsilon
n}\min\\{W\left(\frac{0.45\varepsilon}{\delta}\right),0.6\\}$ where $W(x)$ is
the Lambert W function so $W(x)e^{W(x)}=x$. Given a distribution $P$,

$\mathcal{R}_{\mathcal{N},n,\varepsilon}(P)=\Omega\left(\sum_{x\in[N]}\min\left\\{P(x)(1-P(x)),\sqrt{\frac{P(x)(1-P(x))}{n}}\right\\}+\sum_{x\notin\gamma_{{P}}\left({2\kappa}\right)}P(x)+(|\gamma_{{P}}\left({2\kappa}\right)|-1)\kappa\right)$

Theorem 5.5 follows immediately from Theorem 5.6 and Theorem 5.7. The main
tool we will use is a differentially private version of Assouad’s method. This
gives us a method for lower bounding the error by constructing nets of
distributions that are pairwise far in the relevant metric of interest, which
for us in the TV distance. The following is a slight variant on the
differentially private variant of Assouad’s lemma given in [ASZ21]. Rather
than building a set of distributions indexed by a hypercube, we will build a
set of distributions over a product of hypercubes. Since this is an extension
of the version that appears in [ASZ21], we include a proof in Appendix C for
completeness.

######  Lemma 5.8.

[A extension of $(\varepsilon,\delta)$-DP Assouad’s method [ASZ21]] Let
$k_{0},k_{1},\cdots$ be a sequence of natural numbers such that
$\sum_{s}k_{s}<\infty$, $\varepsilon>0$ and $\delta\in[0,1]$. Given a family
of distributions $\mathcal{P}\subset\Delta(\mathcal{X})$ on a space
$\mathcal{X}$, a parameter $\theta:\mathcal{P}\to\mathcal{M}$ where
$\mathcal{M}$ is a metric space with metric $d$, suppose that there exists a
set $\mathcal{V}\subset\mathcal{P}$ of distributions indexed by the product of
hypercubes $\mathcal{E}_{k_{0}}\times\mathcal{E}_{k_{1}}\times\cdots$ where
$\mathcal{E}_{k}:=\\{\pm 1\\}^{k}$ such that for a sequence
$\tau_{0},\tau_{1},\cdots$,

$\forall(u^{0},u^{1},\cdots),(v^{0},v^{1},\cdots)\in\mathcal{E}_{k_{0}}\times\mathcal{E}_{k_{1}}\times\cdots,\;\;\;d(\theta(p_{u}),\theta(p_{v}))\geq
2\sum_{s}\tau_{s}\sum_{j=1}^{k_{s}}\chi_{u^{s}_{j}\neq v_{j}^{s}}.$ (5)

For each coordinate $s\in\mathbb{N}$, $j\in[k_{s}]$, consider the mixture
distributions obtained by averaging over all distributions with a fixed value
at the $(s,j)$th coordinate:

$p_{+(s,j)}=\frac{2}{|\mathcal{E}_{k_{0}}\times\mathcal{E}_{k_{1}}\times\cdots|}\sum_{u\in\mathcal{E}_{k_{0}}\times\mathcal{E}_{k_{1}}\times\cdots:u^{s}_{j}=+1}p_{u},\;\;p_{-(s,j)}=\frac{2}{|\mathcal{E}_{k_{0}}\times\mathcal{E}_{k_{1}}\times\cdots|}\sum_{u\in\mathcal{E}_{k_{0}}\times\mathcal{E}_{k_{1}}\times\cdots:u^{s}_{j}=-1}p_{u},$

and let $\phi_{s,j}:\mathcal{X}^{n}\to\\{-1,+1\\}$ be a binary classifier.
Then

$\min_{\mathcal{A}\text{ is
}(\varepsilon,\delta)\text{-DP}}\max_{p\in\mathcal{V}}\mathcal{R}_{\mathcal{A},n}(p)\geq\frac{1}{2}\sum_{s}\tau_{s}\sum_{j=1}^{k_{s}}\min_{\phi_{s,j}\text{
is }(\varepsilon,\delta)\text{-DP}}(\Pr_{X\sim
p_{+(s,j)}^{n}}(\phi_{s,j}(X)\neq 1)+\Pr_{X\sim
p_{-(s,j)}^{n}}(\phi_{s,j}(X)\neq-1)),$

where the min on the LHS is over all $(\varepsilon,\delta)$-DP mechanisms, and
on the right hand side is over all $(\varepsilon,\delta)$-DP binary
classifiers. Moreover, if for all $s\in\mathbb{N}$, $j\in[k_{s}]$, there
exists a coupling $(X,Y)$ between $p_{+(s,j)}^{n}$ and $p_{-(s,j)}^{n}$ with
$\mathbb{E}[d_{Ham}(X,Y)]\leq D_{s}$, then

$\min_{\mathcal{A}\text{ is
}(\varepsilon,\delta)\text{-DP}}\max_{p\in\mathcal{V}}\mathcal{R}_{\mathcal{A},n}(p)\geq\sum_{s}\frac{k_{s}\tau_{s}}{2}(0.9e^{-10\varepsilon
D_{s}}-10D_{s}\delta)$

Note that an upper bound on $TV(P_{i},P_{j})\leq\gamma$ implies there exists a
coupling $(X,Y)$ between $P_{i}^{n}$ and $P_{j}^{n}$ such that
$\mathbb{E}[d_{Ham}(X,Y)]\leq n\gamma$.

We will separately prove that each of the three terms in Theorem 5.7 belong in
the lower bound. Each proof will follow the same underlying structure. Given a
distribution $P$, the main technical step is carefully designing a family of
distributions in $\mathcal{N}(P)$ that satisfy the conditions of Lemma 5.8.
Lemma 5.9 and Lemma 5.10 give lower bounds on the noise due to privacy. Lemma
5.11 gives lower bounds based on the error due to sampling.

Let

$\kappa=\frac{1}{10\varepsilon}\min\\{W\left(\frac{0.45\varepsilon}{\delta}\right),0.6\\},$

where $W(x)\approx\ln x-(1-o(1))\ln\ln x$ is the Lambert W function satisfying
$W(x)e^{W(x)}=x$. In both lemma proofs we will use the inequality that if
$D\leq\kappa$, then

$\displaystyle 0.9e^{-10\varepsilon D}-10D\delta$ $\displaystyle\geq
e^{-10\varepsilon
D}\left(0.9-W\left(\frac{0.45\varepsilon}{\delta}\right)e^{W\left(\frac{0.45\varepsilon}{\delta}\right)}\frac{\delta}{\varepsilon}\right)=e^{-10\varepsilon
D}\left(0.9-\frac{0.45\varepsilon}{\delta}\frac{\delta}{\varepsilon}\right)\geq
e^{-10\varepsilon D}0.45\geq 0.2$ (6)

######  Lemma 5.9.

Given a distribution $P$, $\varepsilon>0$, $\delta\in[0,1]$ and
$n\in\mathbb{N}$,

$\mathcal{R}_{\mathcal{N},n,\varepsilon}(P)\geq
0.1\left(\left|\gamma_{{P}}\left({\frac{2\kappa}{n}}\right)\right|-1\right)\frac{\kappa}{n}.$

###### Proof.

Let $L=|\gamma_{{P}}\left({2\kappa/n}\right)|$ be the number of active nodes.
If $L=1$ then the RHS is 0 and so we are done. Otherwise, assume $L>1$ and let
$k={\lfloor{L/2}\rfloor}\geq 1$. Using the notation from Lemma 5.8, let
$k_{0}=k$ and $k_{s}=0$ for all $s>0$. We will drop the reference to $s$ in
the notation since only $s=0$ is significant.

Pair up the active nodes to form $k$ pairs of active nodes denoted by
$(a_{1}^{+},a_{1}^{-}),\cdots,(a_{k}^{+},a_{k}^{-})$. Given
$u\in\mathcal{E}_{k}$, define the distribution $p_{u}$ as follows: for all
$a_{j}^{b}\in\gamma_{{P}}\left({2\kappa/n}\right)$,
$P_{u}(a_{j}^{+})=P(a_{j}^{+})+(\kappa/n)$ and
$P_{u}(a_{j}^{-})=P(a_{j}^{-})-(\kappa/n)$ if $u_{j}=+1$ and
$P_{u}(a_{j}^{+})=P(a_{j}^{+})-(\kappa/n)$ and
$P_{u}(a_{j}^{-})=P(a_{j}^{-})+(\kappa/n)$ if $u_{j}=-1$. For all other $x$,
$P_{u}(x)=P(x)$. It is immediate that for all $u$, $P_{u}\in\mathcal{N}(P)$.
For any pair $u,v$, $d(\theta(p_{u}),\theta(p_{v}))=\text{\rm
TV}(p_{u},p_{v})=d_{Ham}(u,v)(\kappa/n)$, so that Equation 5 is satisfied with
$\tau=\frac{1}{2}(\kappa/n)$. Further, given $j\in[k]$, $p_{+j}$ and $p_{-j}$
only differ on the probability of $a_{j}^{+}$ and $a_{j}^{-}$, so
$D/n=\max_{j}\text{\rm TV}(p_{+j},p_{-j})=\kappa/n$ and by Equation 6,
$0.9e^{-10\varepsilon D}-10D\delta\geq 0.2.$ Noting that
$k\geq(1/2)(\gamma_{{P}}\left({2\kappa/n}\right)-1)$ completes the proof. ∎

######  Lemma 5.10.

For all $\varepsilon>0$, $\delta\in[0,1]$, $n\in\mathbb{N}$ and distributions
$P$ on $[N]$, if $\kappa<n/2$, then

$\mathcal{R}_{\mathcal{N},n,\varepsilon}(P)\geq\Omega\left(\sum_{x\notin\gamma_{{P}}\left({2\kappa}\right)}P(x)\right).$

Since $\kappa\leq\frac{1}{\varepsilon n}$, the condition that $\kappa<n/2$ is
a mild condition. For example, it is satisfied whenever $\varepsilon>2/n^{2}$.

Similar to the proof of Lemma 5.9, we are going to pair up the coordinates and
move mass between the coordinates to create the distributions indexed by the
product of hypercubes. Since we want all the distributions we create to be in
$\mathcal{N}(P)$, we will divide the space into scales such that all elements
in the same scale have approximately the same probability of occurring. We’ll
then move mass within these scales. For $s\in\mathbb{N}$, let
$\mathcal{S}_{s}=\\{x\in[N]\;|\;P(x)\in(2^{-s-1},2^{-s}]\\}$.

###### Proof.

Given $s\in\mathbb{N}$, let
$\mathcal{S}_{s}^{\prime}=\mathcal{S}_{s}\cap\\{x\;|\;P(x)\leq 2\kappa/n\\}$
and $d_{s}=|\mathcal{S}_{s}^{\prime}|$.

Let us first consider the case that there exists a scale $s^{*}$ with
$d_{s^{*}}=1$ and
$P(x^{*})\geq\frac{1}{8}\sum_{x\notin\gamma_{{P}}\left({2\kappa}\right)}P(x)$
where $x^{*}$ is the element in $\mathcal{S}_{s^{*}}^{\prime}$. Define
$P^{\prime}$ by $P^{\prime}(x^{*})=(1/2)P(x^{*})$ and for all $x\neq x^{*}$,
$P^{\prime}(x)=\frac{1-(1/2)P(x^{*})}{1-P(x^{*})}P(x)$. Since
$(1/2)P(x^{*})\leq 2\kappa/n\leq 1/2$, $P^{\prime}\in\mathcal{N}(P)$. In this
case we will use Lemma 5.8 with $k_{0}=1$, $k_{s}=0$ otherwise, and
$\mathcal{E}_{k_{0}}$ corresponds to the set of distributions
$\\{P,P^{\prime}\\}$. Then noting that $\text{\rm
TV}(P,P^{\prime})=(1/2)P(x^{*})$ and using eqn (6) we have that
$\tau=\frac{1}{4}P(x^{*})$ and $D=(1/2)P(x^{*})n\leq\kappa$ so that
$\mathcal{R}_{\mathcal{N},n,\varepsilon}(P)\geq(1/8)P(x^{*})(0.2)=\Omega\left(\sum_{x\notin\gamma_{{P}}\left({2\kappa}\right)}P(x)\right)$
and we are done.

Next suppose that for all scales $s$ such that $d_{s}=1$ we have
$P(x^{*})<\frac{1}{8}\sum_{x\notin\gamma_{{P}}\left({2\kappa}\right)}P(x)$.
Let $s^{*}$ be the smallest $s$ such that $d_{s}=1$. Since the scales
$2^{-s-1}$ are geometrically decreasing,

$\sum_{s:d_{s}=1}\sum_{x\in\mathcal{S}_{s}\cap\\{x|P(x)\leq
2\kappa/n\\}}P(x)\leq 2\sum_{s:d_{s}=1}2^{-s-1}\leq 4\cdot
2^{-s^{*}-1}\leq\frac{1}{2}\sum_{x\notin\gamma_{{P}}\left({2\kappa}\right)}P(x).$

It follows that
$\sum_{s:d_{s}>1}\sum_{x\in\mathcal{S}_{s}\cap\\{x|P(x)<2\kappa/n\\}}P(x)\geq(1/2)\sum_{x\notin\gamma_{{P}}\left({2\kappa}\right)}P(x)$.
Further,

$\sum_{x\notin\gamma_{{P}}\left({2\kappa}\right)}P(x)\leq
2\sum_{s:d_{s}>1}\sum_{x\in\mathcal{S}_{s}\cap\\{x|P(x)<2\kappa/n\\}}P(x)\leq
4\sum_{s:d_{s}>1}2^{-s-1}d_{s}\leq
16\sum_{s:d_{s}>1}2^{-s-1}{\lfloor{d_{s}/2}\rfloor}.$

Thus we can (up to constants) ignore scales such that $d_{s}\leq 1$ and assume
that $d_{s}$ is even for all scales.

Let $k_{s}={\lfloor{d_{s}/2}\rfloor}$. Now, within each scale
$\mathcal{S}_{s}$, pair the elements to form $k_{s}$ distinct pairs
$(a_{s,j}^{+},a_{s,j}^{-})$. Given
$(u^{0},u^{1},\cdots)\in\mathcal{E}_{k_{0}}\times\mathcal{E}_{k_{1}}\times\cdots$,
define $p_{u}$ by $p_{u}(a_{s,j}^{+})=p(a_{s,j})+2^{-s-2}$ and
$p_{u}(a_{s,j}^{-})=p(a_{s,j})-2^{-s-2}$ if $u^{s}_{j}=+1$ and
$p_{u}(a_{s,j}^{+})=p(a_{s,j})-2^{-s-2}$ and
$p_{u}(a_{s,j}^{-})=p(a_{s,j})+2^{-s-2}$ if $u^{s}_{j}=-1$. For all other
elements, $p_{u}(x)=p(x)$. Then, it is easy to see that for all $u$,
$p_{u}\in\mathcal{N}(P)$. Further, using notation from Lemma 5.8, Equation 5
is satisfied with $\tau_{s}=\frac{1}{2}2^{-s-2}$ since
$d(\theta(p_{u}),\theta(p_{v}))=\text{\rm
TV}(p_{u},p_{v})=2\sum_{s}\tau_{s}d_{Ham}(u^{s},v^{s})$ and
$D_{s}/n=\max_{j}\text{\rm TV}(p_{+(s,j)},p_{-(s,j)})=2^{-s-2}$, which is less
than $\kappa$ whenever $k_{s}>0$. By eqn (6), $0.9e^{-10\varepsilon
D_{s}}-10D_{s}\delta\geq 0.2$ whenever $k_{s}>0$ so by Lemma 5.8 we have

$\mathcal{R}_{\mathcal{N},n,\varepsilon}(P)\geq\sum_{s}\frac{1}{2}2^{-s-2}k_{s}(0.2)=\frac{0.2}{4}\sum_{s:d_{s}>1}2^{-s-1}{\lfloor{d_{s}/2}\rfloor}\geq\frac{0.2}{4\times
16}\sum_{x\notin\gamma_{{P}}\left({2\kappa}\right)}P(x)$

which completes the proof. ∎

Next we lower bound the statistical term.

######  Lemma 5.11.

For all $n\in\mathbb{N}$, $\varepsilon>0$, $\delta\in[0,1]$ and distributions
$P$, if $n\geq 2$ and $\varepsilon>2/n$, then

$\mathcal{R}_{\mathcal{N},n,\varepsilon}(P)\geq\mathcal{R}_{\mathcal{N},n}(P)\geq\Omega\left(\sum_{x\in[N]}\min\left\\{P(x)(1-P(x)),\sqrt{\frac{P(x)(1-P(x)}{n}}\right\\}\right).$

To streamline the notation, we will use $L(x)$ to denote
$\min\left\\{x(1-x),\sqrt{\frac{x(1-x)}{n}}\right\\}$. In order to prove Lemma
5.11, we will need the following standard result from the statistics
literature which allows us to lower bound the performance of any simple
classifier distinguishing two distributions $P$ and $Q$ by the KL divergence
between $P$ and $Q$. We give a specific result for distinguishing Bernoulli
random variables since we’ll use this in the proof of Lemma 5.11.

######  Lemma 5.12.

Given any pair of distributions $P$ and $Q$ on the same domain,

$\min_{\phi}\left(\Pr_{X\sim P^{n}}(\phi(X)=1)+\Pr_{X\sim
Q^{n}}(\phi(X)=-1)\right)\geq\frac{1}{2}(1-\sqrt{n\text{\rm KL}(P,Q)}),$

where the minimum is over all binary classifiers. In particular, if
$P=\texttt{Bernoulli}(p-\alpha)$ and $Q=\texttt{Bernoulli}(p+\alpha)$ where
$0\leq\alpha\leq\frac{1}{2}L(p)$ then

$\min_{\phi}\left(\Pr_{X\sim P^{n}}(\phi(X)=1)+\Pr_{X\sim
Q^{n}}(\phi(X)=-1)\right)\geq 1/4,$

where again the minimum is over all binary classifiers.

The proof of Lemma 5.12 can be found in Appendix C

###### Proof of Lemma 5.11.

As in the proof of Lemma 5.10, first suppose there exists a scale $s^{*}$ with
$d_{s^{*}}=1$ and there exists $x^{*}\in\mathcal{S}_{s^{*}}$ such that

$\frac{1}{2}L(P(x^{*}))\geq\frac{1}{60}\sum_{x\in[N]}L(P(x)).$

Then define a distribution $P^{\prime}$ by
$P^{\prime}(x^{*})=P(x^{*})-\frac{1}{2}L(P(x^{*}))$ and for all $x\neq x^{*}$,
$P^{\prime}(x)=\frac{1-P(x^{*})+\frac{1}{2}L(P(x^{*}))}{1-P(x^{*})}P(x)$. Then
$P^{\prime}\in\mathcal{N}(P)$ since
$\frac{1}{2}L(P(x^{*}))<\frac{1}{2}\min\\{P(x^{*}),(1-P(x^{*}))\\}$. Then we
will use Lemma 5.8 with $k_{0}=1$ and $k_{s}=0$ for $s>0$, and
$\mathcal{E}_{k_{0}}$ corresponds to $\\{P,P^{\prime}\\}$. Now,

$\displaystyle\text{\rm KL}(P^{\prime},P)$
$\displaystyle=(P(x^{*})-\frac{1}{2}L(P(x^{*})))\ln\frac{P(x^{*})-\frac{1}{2}L(P(x^{*}))}{P(x^{*})}+(1-P(x^{*})+\frac{1}{2}L(P(x^{*})))\ln\frac{1-P(x^{*})+\frac{1}{2}L(P(x^{*}))}{1-P(x^{*})}$
$\displaystyle\leq\frac{1}{4n}$

(for more detail on the proof of this inequality see the proof of Lemma 5.12)
so

$\min_{\phi}\left(\Pr_{X\sim
P^{n}}(\phi(X)=1)+\Pr_{X\sim{P^{\prime}}^{n}}(\phi(X)=-1)\right)\geq\frac{1}{2}(1-\sqrt{n\text{\rm
KL}(P,P^{\prime})})\geq 1/4$

and $\tau_{0}=\text{\rm TV}(P,P^{\prime})=\frac{1}{2}L(P(x^{*}))$. Thus by
Lemma 5.8,

$\mathcal{R}_{\mathcal{N},n,\varepsilon}(P)\geq\mathcal{R}_{\mathcal{N},n}(P)\geq\frac{1}{2}L(P(x^{*}))\frac{1}{4}\geq\frac{1}{480}\sum_{x\in[N]}\frac{1}{2}L(P(x)),$

and we are done.

On the other hand, suppose that for all scales $s$ such that $d_{s}=1$ we have

$L(P(x_{s}))~{}\leq~{}\frac{1}{30}\sum_{x\in[N]}L(P(x)),$

where $\mathcal{S}_{s}=\\{x_{s}\\}$. As in the proof of Lemma 5.10, we will
argue that we can ignore any singleton scales, and assume that $d_{s}$ is even
for all scales. Let $s^{*}=\min\\{s>0\;|\;d_{s}=1\\}$ so

$\displaystyle\sum_{s:d_{s}=1}L(P(x_{s}))$
$\displaystyle\leq\chi_{d_{0}=1}L(P(x_{0}))+\sum_{s>0:d_{s}=1}\min\left\\{2^{-s},\sqrt{\frac{2^{-s}}{n}}\right\\}$
$\displaystyle\leq\chi_{d_{0}=1}L(P(x_{0}))+(2+\sqrt{2})\min\left\\{2^{-s^{*}},\sqrt{\frac{2^{-s^{*}}}{n}}\right\\}$
$\displaystyle\leq\chi_{d_{0}=1}L(P(x_{0}))+2(2+\sqrt{2})\min\left\\{P(x_{s^{*}}),\sqrt{\frac{P(x_{s^{*}})}{n}}\right\\}$
$\displaystyle\leq\chi_{d_{0}=1}L(P(x_{0}))+4(2+\sqrt{2})L(P(x_{s^{*}}))$
$\displaystyle\leq\frac{1+4(2+\sqrt{2})}{30}\sum_{x\in[N]}L(P(x)).$

Therefore,
$\sum_{s:d_{s}>1}\sum_{x\in\mathcal{S}_{s}}L(P(x))\geq(1/2)\sum_{x\in[N]}L(P(x))$
and so

$\displaystyle\sum_{s}L(2^{-s-1}){\lfloor{d_{s}/2}\rfloor}$
$\displaystyle\geq\sum_{s:d_{s}>1}L(2^{-s-1})\frac{1}{3}d_{s}$
$\displaystyle\geq\sum_{s:d_{s}>1}\sum_{x\in\mathcal{S}_{s}}L(2^{-s-1})\frac{1}{3}$
$\displaystyle\geq\frac{1}{3\sqrt{2}}\sum_{x\in[N]}L(P(x))$ (7)

where the first inequality follows from
${\lfloor{d_{s}/2}\rfloor}\geq(1/3)d_{s}$ whenever $d_{s}>1$, and the second
follows because $2^{-s-1}\leq P(x)\leq 1/2$ for all $x\in\mathcal{S}_{s}$ such
that $d_{s}>1$.

Assume that $d_{s}$ is even for all $s$. Within each scale $\mathcal{S}_{s}$,
pair the elements to form $k_{s}=d_{s}/2$ distinct pairs
$(a^{+}_{s,j},a^{-}_{s,j})$ per scale. For all $s\in\mathbb{N}$, let
$\alpha_{s}=\frac{1}{2}L(2^{-s-1})$, and note that for all
$x\in\mathcal{S}_{s}$ and $s>0$, $\alpha_{s}\leq\frac{1}{2}L(P(x))$. Given
$(u^{0},u^{1},\cdots)\in\mathcal{E}_{k_{0}}\times\mathcal{E}_{k_{1}}\times\cdots$,
define $p_{u}$ by $p_{u}(a_{s,j}^{+})=p(a_{s,j}^{+})+\alpha_{s}$ and
$p_{u}(a_{s,j}^{-})=p(a_{s,j}^{-})-\alpha_{s}$ if $u^{s}_{j}=+1$ and
$p_{u}(a_{s,j}^{+})=p(a_{s,j}^{+})-\alpha_{s}$ and
$p_{u}(a_{s,j}^{-})=p(a_{s,j}^{-})+\alpha_{s}$ if $u^{s}_{j}=-1$. For all
other elements, $p_{u}(x)=p(x)$. Then, for all $u$, $p_{u}\in\mathcal{N}(P)$.
Further, using notation from Lemma 5.8, we have $\tau_{s}=\alpha_{s}$. Also,
for any $(s,j)$, $p_{+(s,j)}$ and $p_{-(s,j)}$ only differ on $a^{+}_{s,j}$
and $a^{-}_{s,j}$ where $p_{+(s,j)}(a^{+}_{s,j})=P(a^{+}_{s,j})+\alpha_{s}$
and $p_{-(s,j)}(a^{+}_{s,j})=P(a^{+}_{s,j})-\alpha_{s}$. Therefore, by Lemma
5.12, and the post-processing inequality,

$\min_{\phi}\left(\Pr_{X\sim{p_{+(s,j}}^{n}}(\phi(X)=1)+\Pr_{X\sim{p_{-(s,j}}^{n}}(\phi(X)=-1)\right)\geq
1/4.$

Lemma 5.8 then implies the result. ∎

Theorem 5.7 follows immediately from Lemma 5.9, Lemma 5.10 and Lemma 5.11.

###  5.3 An $\varepsilon$-DP Distribution Estimation Algorithm

Now, let us return to HSTs and designing an estimation algorithm that achieves
the target estimation rate, up to logarithmic factors. As in the one-
dimensional setting, we want to restrict to only privately estimating the
density at a small number ($\approx\varepsilon n$) of points. While we could
try to mimic the one-dimensional solution by privately estimating a solution
to the $\varepsilon n$-median problem, it’s not clear how to prove that such
an approach is instance-optimal. It turns out that a simpler solution more
amenable to analysis will suffice. Our algorithm has two stages; first we
attempt to find the set of $\frac{\log(1/\delta)}{\varepsilon n}$-active
nodes, then we estimate the weight of these active nodes. Since these nodes
have weight greater than $\frac{\log(1/\delta)}{\varepsilon n}$, we can
privately estimate them to within constant multiplicative error.

Let $\mathcal{X}$ be the underlying metric space so $P\in\Delta(\mathcal{X})$.
For any set $S$ of nodes and a function $F$ defined on the nodes, define the
function $F|_{S}$ as $F|_{S}(\nu)=F(\nu)$ if $\nu\in S$ and $F|_{S}(\nu)=0$
otherwise. Given two functions $F$ and $G$ defined on the nodes, we define

$\mathfrak{W}(F,G)=\sum_{\nu}r_{\nu}|F(\nu)-G(\nu)|,$

where $r_{\nu}$ is the length of the edge connecting $\nu$ to its parent, and
the sum is over all nodes in the tree. So by Lemma 5.4,
$\mathcal{W}(P,Q)=\mathfrak{W}(\mathfrak{G}_{P},\mathfrak{G}_{Q})$. Note that
$\mathfrak{W}$ satisfies the triangle inequality.

Algorithm 1 PrivDensityEstTree

1:Input: $D\in\mathcal{X}^{n},\varepsilon$

2:$\widehat{\mathfrak{G}_{P}}=\texttt{EmpDist}(D)$ $\triangleright$ Compute
empirical distribution.

3:$\hat{\gamma}_{\varepsilon}=\texttt{LocateActiveNodes}(\widehat{\mathfrak{G}_{P}};\varepsilon)$
$\triangleright$ Privately approximate set of active nodes.

4:Define $\widetilde{\mathfrak{G}_{\hat{P}_{n},\hat{\gamma}_{\varepsilon}}}$
by
$\widetilde{\mathfrak{G}_{\hat{P}_{n},\hat{\gamma}_{\varepsilon}}}(x)=\begin{cases}0&\text{if
}x\notin\hat{\gamma}_{\varepsilon}\\\
\widehat{\mathfrak{G}_{P}}(x)+\mathsf{Lap}(\frac{1}{\varepsilon
n}))&\text{otherwise.}\end{cases}$$\triangleright$ Approximate densities.

5:$\hat{P}_{n,\varepsilon}=\texttt{Projection}(\widetilde{\mathfrak{G}_{\hat{P}_{n},\hat{\gamma}_{\varepsilon}}})$
$\triangleright$ Project noisy densities onto space of distributions.

6:return $\hat{P}_{n,\varepsilon}$

A high-level outline of the proposed algorithm is given in Algorithm 1. Now,
we state the main theorem of this section.

######  Theorem 5.13.

Given any $\varepsilon>0$, PrivDensityEstTree is
$(0pt+1)\varepsilon$-DP. Given a distribution $P$, with probability
$1-(0pt\log n+40pt\varepsilon n)\beta$,

$\displaystyle\mathcal{W}(P,\hat{P}_{\varepsilon})=O\Bigg{(}\sum_{\ell\in[0pt]}\sum_{x\in[N_{\ell}]}\min\left\\{P_{\ell}(x),1-P_{\ell}(x)\sqrt{\frac{P_{\ell}(x)\log(n/\beta)}{n}}\right\\}$
$\displaystyle\hskip
72.26999pt+\sum_{\nu\notin\gamma_{{P}}\left({\max\\{\frac{2}{\varepsilon
n}+2\frac{\log(2/\beta)}{\varepsilon n},\frac{192\log(n/\beta)}{\varepsilon
n}\\}}\right)}\mathfrak{G}_{P}(\nu)+\frac{|\gamma_{{P_{\ell}}}\left({\frac{1}{2\varepsilon
n}}\right)-1|\log(1/\beta)}{\varepsilon n}\Bigg{)}$

This bound has the same three terms as our lower bound on
$\mathcal{R}_{\mathcal{N},n,\varepsilon}$ in Theorem 5.5 corresponding again
to the empirical error (the error inherent even in the absence of a privacy
requirement), the error from the private algorithm not being able to estimate
the probability of events that occur with probability less than
$\approx\log(1/\delta)/\varepsilon n$, and the error due to the noise added to
the active nodes. The maximum over the levels that appeared in the lower bound
is replaced with a sum over the levels in the upper bound, so, up to
logarithmic factors, the upper bound is within a factor of $0pt$ of the lower
bound. Since we can not hope to locate the set of $\log(1/\delta)/(\varepsilon
n)$-active nodes exactly with a private algorithm, we find a set
$\hat{\gamma}_{n}$ that is guaranteed to satisfy

$\gamma_{{P}}\left({\max\left\\{\frac{2}{\varepsilon
n}+2\frac{\log(2/\beta)}{n},\frac{192\log(n/\beta)}{n}\right\\}}\right)\subset\hat{\gamma}_{\varepsilon}\subset\gamma_{{P}}\left({\frac{1}{2\varepsilon
n}}\right).$

Note that $\max\left\\{\frac{2}{\varepsilon
n}+2\frac{\log(2/\beta)}{n},\frac{192\log(n/\beta)}{n}\right\\}\leq\frac{C\log(n/\beta)}{\varepsilon
n}$ so the error introduced here by not estimating
$\gamma_{{P}}\left({\frac{1}{\varepsilon n}}\right)$ perfectly is at most a
logarithmic multiplicative factor.

Algorithm 2 EmpDist

1:Input: $D\in\mathcal{X}^{n},A$

2:Let $\hat{P}_{n}$ be the empirical distribution.

3:for all node $\nu$ do

4:
$\widehat{\mathfrak{G}_{P}}(\nu)=\begin{cases}0&\mathfrak{G}_{\hat{P}_{n}}(\nu)<\frac{\sqrt{\log(n/\beta)}}{n}\\\
1&\mathfrak{G}_{\hat{P}_{n}}(\nu)>1-\frac{\sqrt{\log n/\beta)}}{n}\\\
\mathfrak{G}_{\hat{P}_{n}}(\nu)&\text{otherwise}\end{cases}$

The first step of our algorithm is to estimate the empirical distribution. We
use a truncated version of the standard empirical distribution. This allows us
to achieve an error rate of $~{}\min\\{P(x),\sqrt{P(x)/n}\\}$ even when $P(x)$
is small.

The proof of the following lemma is contained in Appendix C.

######  Lemma 5.14.

For any distribution $P$, if $\log(n/\beta)>1$ then with probability
$1-30pt\beta$,

$\mathfrak{W}(\widehat{\mathfrak{G}_{P}},\mathfrak{G}_{P})\leq\sum_{\ell\in[0pt]}\sum_{x\in[N_{\ell}]}\min\left\\{P_{\ell}(x)(1-P_{\ell}(x)),4\sqrt{3\frac{P_{\ell}(x)(1-P_{\ell}(x))\log(n/\beta)}{n}}\right\\}$

Algorithm 3 LocateActiveNodes

1:Input: $\widehat{\mathfrak{G}_{P}},\varepsilon$

2:Let $\ell=0$ and $\hat{\gamma}_{\varepsilon,0}=\\{\nu\\}$ where $\nu$ is the
root node.

3:while $\hat{\gamma}_{\varepsilon,\ell}\neq\varnothing$ and $\ell<0pt$ do

4: $\hat{\gamma}_{\varepsilon,\ell+1}=\varnothing$

5: for all $\nu\in\hat{\gamma}_{\varepsilon,\ell}$ do

6: for all children $\nu^{\prime}$ of $\nu$ do

7: if
$\widehat{\mathfrak{G}_{P}}(\nu^{\prime})+\mathsf{Lap}(\frac{1}{\varepsilon
n})>2\kappa+\frac{\log(2/\beta)}{\varepsilon n}$ then

8:
$\hat{\gamma}_{\varepsilon,\ell+1}=\hat{\gamma}_{\varepsilon,\ell+1}+\\{\nu^{\prime}\\}$

9: $\ell=\ell+1$

10:return $\cup\hat{\gamma}_{\varepsilon,\ell}$

The goal of Algorithm 3 is to estimate the set of $1/(\varepsilon n)$-active
nodes.

The next lemma allows us to bound how close to the goal we get. The proof is
contained in Appendix C.

######  Lemma 5.15.

Let $\hat{\gamma}_{\varepsilon}$ be the set of active nodes found in Algorithm
1. Then with probability $1-0pt(\log n+4\varepsilon n)\beta$,

$\gamma_{{P}}\left({\max\left\\{\frac{2}{\varepsilon
n}+4\frac{\log(2/\beta)}{\varepsilon
n},\frac{192\log(n/\beta)}{n}\right\\}}\right)\subset\hat{\gamma}_{\varepsilon}\subset\gamma_{{P}}\left({\frac{1}{2\varepsilon
n}}\right).$

We also prove the following lemma relating the error due to estimating the
active nodes to a quantity depending on the true active nodes.

######  Lemma 5.16.

If $\gamma_{{P}}\left({\max\\{\frac{2}{\varepsilon
n}+4\frac{\log(2/\beta)}{\varepsilon
n},\frac{192\log(n/\beta)}{n}\\}}\right)\subset\hat{\gamma}_{\varepsilon}$
then

$\mathfrak{W}(\widehat{\mathfrak{G}_{P}},\widehat{\mathfrak{G}_{P}}|_{\hat{\gamma}_{\varepsilon}})\leq\mathfrak{W}(\mathfrak{G}_{P},\widehat{\mathfrak{G}_{P}})+\mathfrak{W}(\mathfrak{G}_{P},\mathfrak{G}_{P}|_{\gamma_{{P}}\left({\max\\{\frac{2}{\varepsilon
n}+4\frac{\log(2/\beta)}{\varepsilon
n},\frac{192\log(n/\beta)}{n}\\}}\right)})$

Algorithm 4 Projection

1:Input: $\mathfrak{G}$, a real-valued function on the nodes of the HST such
that $\mathfrak{G}(\nu_{0})=1$ where $\nu_{0}$ is the root node.

2:$\bar{\mathfrak{G}}=\mathfrak{G}$

3:for $\ell=0:0pt-1$ do

4: for all nodes $\nu$ at level $\ell$ do

5: Let $A_{\nu}=\sum\mathfrak{G}(\nu^{\prime})$ where the sum is over the
children of $\nu$.

6: Let $d_{\nu}$ be the number of children of $\nu$

7: if $A_{\nu}=0$ then

8: for all children $\nu^{\prime}$ of $\nu$ do

9: $\bar{\mathfrak{G}}(\nu^{\prime})=\frac{1}{d_{\nu}}\bar{\mathfrak{G}}(\nu)$

10: else

11: for all children $\nu^{\prime}$ of $\nu$ do

12:
$\bar{\mathfrak{G}}(\nu^{\prime})=\frac{\bar{\mathfrak{G}}(\nu)}{A_{\nu}}\mathfrak{G}(\nu^{\prime})$

13:return $\bar{\mathfrak{G}}$

The key component of this proof is that any discrepancy between the weight of
the nodes on $P$ and that assigned by $\widehat{\mathfrak{G}_{P}}$ was already
paid for in $\mathcal{W}(P,\widehat{\mathfrak{G}_{P}})$. The final step in
Algorithm 1 is to project the noisy function
$\widetilde{\mathfrak{G}_{\hat{P_{n}},\hat{\gamma_{\varepsilon}}}}$ into the
space of distributions on the underlying metric space. We’d like to do this in
a way that preserves, up to a constant, the $\mathfrak{W}$ distance between
$P$ and $\widetilde{\mathfrak{G}_{\hat{P_{n}},\hat{\gamma_{\varepsilon}}}}$.
We will do this iteratively starting from the root node, by ensuring that the
sum of each node’s children add up to it’s assigned value. Since we know the
root node has value 1, this results in a valid distribution. We start from the
top of the tree since errors in higher nodes of the contribute more to the
Wasserstein distance. While errors in higher nodes of the tree propagate can
propagate to lower levels, the predominant influence on the overall error is
retained at the top level due to the geometric nature of the edge weights.

######  Lemma 5.17.

For any real-valued function $\mathfrak{G}$ on the nodes of the HST such that
$\mathfrak{G}(\nu_{0})=1$ where $\nu_{0}$ is the root node and given any
distribution $P$,

$\mathcal{W}(P,\texttt{Projection}(\mathfrak{G}))\leq
4\mathfrak{W}(\mathfrak{G}_{P},\mathfrak{G}).$

Combining the above lemmas appropriately gives the proof of Theorem 5.13 (see
Appendix C).

###### Proof of Theorem 5.13.

The privacy follows from the fact that each user contributes to at most $0pt$
queries in LocateActiveNodes and at most one coordinate in the computation of
$\widetilde{\mathfrak{G}_{\hat{P_{n}},\hat{\gamma}_{n}}}$ in line 4 in
PrivDensityEstTree.

For the utility, we will consider each level $\ell$ individually. First
suppose that $|\gamma_{{P_{\ell}}}\left({1/(2\varepsilon n)}\right)|>1$.

$\displaystyle\mathcal{W}(P_{\ell},(\hat{P}_{\varepsilon})_{\ell})$
$\displaystyle\leq
2\mathfrak{W}((\mathfrak{G}_{P})_{\ell},(\widetilde{\mathfrak{G}_{\hat{P}_{n},\hat{\gamma}_{\varepsilon}}})_{\ell})$
$\displaystyle\leq
2\left(\mathfrak{W}((\mathfrak{G}_{P})_{\ell},(\widehat{\mathfrak{G}_{P}})_{\ell})+\mathfrak{W}((\widehat{\mathfrak{G}_{P}})_{\ell},(\widehat{\mathfrak{G}_{P}}|{\hat{\gamma}_{\varepsilon}})_{\ell})+\mathfrak{W}((\widehat{\mathfrak{G}_{P}}|{\hat{\gamma}_{\varepsilon}})_{\ell},(\widetilde{\mathfrak{G}_{\hat{P}_{n},\hat{\gamma}_{\varepsilon}}})_{\ell})\right)$
(8) $\displaystyle\leq
2\left(2\mathfrak{W}((\mathfrak{G}_{P})_{\ell},(\widehat{\mathfrak{G}_{P}})_{\ell})+\mathfrak{W}((\mathfrak{G}_{P})_{\ell},(\mathfrak{G}_{P}|_{\gamma_{{P}}\left({\max\\{\frac{2}{\varepsilon
n}+2\frac{\log(2/\beta)}{n},\frac{\log(n/\beta)}{n}\\}}\right)})_{\ell})+\mathfrak{W}((\widehat{\mathfrak{G}_{P}}|_{\hat{\gamma}_{\varepsilon}})_{\ell},(\widetilde{\mathfrak{G}_{\hat{P}_{n},\hat{\gamma}_{\varepsilon}}})_{\ell})\right)$

where the first inequality follow from Lemma 5.17, the second inequality
follows from the triangle inequality and Lemma 5.4, and the third follows from
Lemma 5.16 and Lemma 5.15. Finally,

$\mathfrak{W}(\widehat{\mathfrak{G}_{P}}|_{\hat{\gamma}_{\varepsilon}})_{\ell},(\widetilde{\mathfrak{G}_{\hat{P}_{n},\hat{\gamma}_{\varepsilon}}})_{\ell})\leq\sum_{\nu\in\gamma_{{P}}\left({\frac{1}{2\varepsilon
n}}\right)}r_{\nu}|\mathsf{Lap}(\frac{1}{\varepsilon
n})|\leq\frac{1}{2}\sum_{\nu\in\gamma_{{P}}\left({\frac{1}{2\varepsilon
n}}\right)}r_{\nu}|\mathsf{Lap}(\frac{1}{\varepsilon n})|$

The final statement then follows from Lemma 5.14 and basic concentration
bounds on the Laplacian distribution.

If $|\gamma_{{P_{\ell}}}\left({1/(2\varepsilon n)}\right)|=1$, then the proof
goes through for all except the final term related to the noise due to
privacy. We consider two cases. Let
$x\in\gamma_{{P_{\ell}}}\left({1/(2\varepsilon n)}\right)$. First suppose that
$P_{\ell}(x)>1-\frac{1}{2\varepsilon n}$ then no node that is in a level above
$x$, but is not a direct ancestor of $x$ is in
$\gamma_{{P_{\ell}}}\left({1/(2\varepsilon n)}\right)$. Therefore, since the
projection algorithm is top-down, $(\hat{P}_{n,\varepsilon})_{\ell}$ will be
concentrated on $x$. Therefore, the error of level $\ell$ is simply
$(1-P(x))$, which can be charged to the first term plus the sum of the weight
of the inactive nodes, which is in the second term. Next, suppose that
$P_{\ell}(x)<1-\frac{1}{2\varepsilon n}$ then sum of the inactive nodes (in
term two) dominates the error due to adding noise to $P(x)$ ∎

##  6 Instance Optimal Density Estimation on $\mathbb{R}$ in Wasserstein
distance

Let us now consider the setting of estimating distributions $P$ on
$\mathcal{X}=\mathbb{R}$. In this setting, the target estimation rate is that
of an algorithm that knows that the distribution is either $P$ or $Q_{P}$ for
a distribution $Q_{P}$ such that $D_{\infty}(P,Q_{P})\leq\ln 2$. This
definition of instance-optimality strengthens that corresponding to the so-
called hardest-one dimensional subproblem [DL91], since this is a harder
estimation rate to achieve. A formal description of the target estimation rate
is given in Section 3.1 and Section 3.3. In Section 6.1, we lower bound this
estimation rate using hypothesis testing techniques. Then, in Section 6.2, we
give an algorithm that up to polylogarithmic factors, uniformly achieves the
lower bound, and hence approximately achieves the instance-optimal estimation
rate. Our instance optimality results apply to all continuous distributions in
a bounded interval with density functions (though it is likely that they apply
more generally). All omitted proofs can be found in Appendix F.

###  6.1 General Lower Bound

To state the main theorem in this section, we will introduce some notation. We
start by defining the restriction of a distribution.

######  Definition 6.1.

For any distribution $P$ over $\mathbb{R}$ with a density function, the
restriction $P|_{u,v}$ of $P$ with respect to $u\leq v\in\mathbb{R}$ is
defined as the distribution with the following CDF function F’:

$F^{\prime}_{P_{u,v}}(t)=\begin{cases}0&\text{$t<u$}\\\ F_{P}(x)&\text{$u\leq
t<v$}\\\ 1&\text{$t\geq v$}\\\ \end{cases}$

If $u=v$, then $F^{\prime}$ is a step function that goes to $1$ at that point
and is $0$ prior to that point.

Also recall the following definition of quantiles.

######  Definition 6.2.

For $0<\alpha\leq 1$, the $\alpha$-quantile of a distribution $P$ over
$\mathbb{R}$ is defined as follows:

$q_{\alpha}(P)=\arg\min_{t}\\{\Pr_{y\sim P}(y\leq t)\geq\alpha\\}.$

When the distribution $P$ is clear from context, we will sometimes abuse
notation and use $q_{\alpha}$ when we mean $q_{\alpha}(P)$. The main theorem
we will prove in this section is the following:

######  Theorem 6.3.

There exists a constant $C$ such that given a continuous distribution $P$ on
$\mathbb{R}$ with bounded expectation and
$\varepsilon\in(0,1],n\in\mathbb{N}$,

$\displaystyle\mathcal{R}_{loc,n,\varepsilon}(P)=\Omega\Bigg{(}\frac{1}{\varepsilon
n}$ $\displaystyle\left(q_{1-\frac{1}{C\varepsilon
n}}-q_{\frac{1}{C\varepsilon
n}}\right)+\mathcal{W}(P,P|_{q_{\frac{1}{C\varepsilon
n}},q_{1-\frac{1}{C\varepsilon n}}})$ $\displaystyle+\frac{1}{\sqrt{\log
n}}\mathbb{E}\left[\mathcal{W}(P|_{q_{\frac{1}{C\varepsilon
n}},q_{1-\frac{1}{C\varepsilon n}}},\hat{P}_{n}|_{q_{\frac{1}{C\varepsilon
n}},q_{1-\frac{1}{C\varepsilon n}}})\right]\Bigg{)},$

where $\hat{P}_{n}$ is the empirical distribution on $n$ samples drawn
independently from $P$.

The same result can be extended to $(\varepsilon,\delta)$-DP algorithms as
well for $\delta=o(\frac{1}{n})$

We discuss each of the terms in turn. Note that the final term is related to
the expected Wasserstein distance between the empirical distribution and the
true distribution. There is now a long line of work characterizing this
quantity in terms of the distribution (See Section 4), but essentially, if the
distribution is more concentrated, this term is smaller. The first term is a
very particular inter-quantile distance that is also much smaller for
concentrated distributions, and can be large for relatively dispersed
distributions. The second term characterizes the length of the tails of the
distribution—longer tails make this Wasserstein distance larger. Overall, this
rate is significantly lower for more concentrated distributions with small
support, and relatively large for more dispersed distributions. We prove this
theorem over the following couple of sections; in Section 6.1.1 we
characterize the cost of private instance optimality, and in Section 6.1.2 we
characterize the cost of achieving instance optimality without privacy (this
non-private characterization is also new to our work, to the best of our
knowledge). Combining the theorems in those sections gives the above result.

####  6.1.1 The Privacy Term

The main theorem we will prove in this section is the following.

######  Theorem 6.4.

Fix $\varepsilon\in(0,1]$, $n\in\mathbb{N}$. For all distributions $P$ over
$\mathbb{R}$ that have a density function and finite expectation, there exists
another distribution $Q^{\prime\prime}$ such that $D_{\infty}(P,Q)\leq 2$,
that is indistinguishable from $P$ given $O(n)$ samples such that for all
$\varepsilon$-DP algorithms $A:\mathbb{R}^{n}\to\Delta(\mathbb{R})$, with
probability at least $0.25$ over the draws ${\bf x}\sim P^{n}$, ${\bf
x}^{\prime}\sim Q^{\prime\prime n}$, the following holds for some constant
$C$.