Instance-Optimal Private Density Estimation in the Wasserstein Distance

Authors: vitaly  feldman
audra  mcmillan
satchit  sivakumar
kunal  talwar
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 01:31:51

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 slightly 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 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 learning in TV distance for discrete distributions.

Main paper: 
 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, ℓ] 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, ℓ] 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 π from P to Q of the quantity E x∼P [d(x, π(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 Appendix C) . 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) = 1(x ̸ = y), it reduces to the commonly used total variation (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 Pn and the target distribution P over R d . Similarly, this question for the case of the discrete metric (d(x, y) = 1(x ̸ = 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 ε-privately learning a discrete distribution over {0, . . . , 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 ℓ 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 Θ(n -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 Θ((εn) -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 (See Appendix F for details on this experiment). 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 Re d lies on a k-dimensional subspace, the worst-case error scaling with Õ((εn) -1 d ) is significantly larger than our algorithm's scaling of Õ((εn) -1 k ). This motivates the problem of viewing this question through the lens of instance optimality. 2 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 M be a class of algorithms of interest (e.g. all (ε, δ)-differentially private algorithms) and cost(•, 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 Pn of n samples from P , we want to learn an estimate A( Pn ) for the distribution. Our measure of performance is the Wasserstein distance W, so cost(A, P ) = E[W(P, A( Pn ))]. We would ideally like to say that an algorithm A is α-instance optimal in a class M if for all instances P , and all A ′ ∈ M , cost(A( Pn )), P ) ≤ α • cost(A ′ ( Pn )), P ).(InstanceOptimality-Ideal)The reader would have noticed that this definition is however impossible to achieve except for trivial classes M. The algorithm A ′ 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 A ′ not just on the given instance, but on a small neighborhood around it. Thus we say (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.that that an algorithm A is α-instance optimal amongst a class M with respect to a neighborhood function N if for all instances P , and all A ′ ∈ M cost(A( Pn )), P ) ≤ α • supP ′ ∈N (P ) cost(A ′ ( P ′ n )), P ′ ).In other words, the benchmark we evaluate against is the cost of the best algorithm for a neighborhood N (P ) that knows this neighborhood. We would like our algorithm A, that is not tailor-made for 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 N (P ) must be carefully defined for this notion to be meaningful; we can always define 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 N (P ). For example, [FLN01, ABC17, VV16, OS15, GKN20] include in 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 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 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 N (P ). We discuss this general definition of instance optimality further in Appendix B. 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 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 ℓ 1 and ℓ 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 Appendix B and Appendix C 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 ∞ (P, Q) is defined as sup x∈X max ln P (x) Q(x) , ln Q(x) P (x) . Our neighborhood map N will have the property that for all P , and for all Q ∈ N (P ), D ∞ (P, Q) ≤ ln 2. This corresponds to the benchmark algorithm A ′ 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 N (P ) is defined to be {P, Q} where Q is a distribution with D ∞ (P, Q) ≤ 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 Appendix B. Informally, local minimality says that if any comparator algorithm does better than A on P , then there is a distribution Q in the neighborhood of P where 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 Appendix B) 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 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 .Preliminaries: First, we define differential privacy. Further discussion on differential privacy can be found in Appendix A. Definition 1.1 (Differential Privacy [DMNS17, DKM + 06]). A randomized algorithm A : X n → Y is (ε, δ)-differentially private if for every pair of datasets x, x ′ ∈ X n that differ in at-most one data entry, and for all events Y ⊆ Y,Pr[A(x) ∈ Y ] ≤ e ε • Pr[A(x ′ ) ∈ Y ] + δ.Given an estimation algorithm A : X n → M, the estimation rate of A for distribution P is: R A,n (P ) = inf t∈R {t : w.p. ≥ 0.75 over x ∼ P n and the randomness of the algorithm, W d (A(x), P ) ≤ t}.(1)We start by stating an informal version of our result in the one-dimensional real case.Theorem 1.2 (Informal 1-dimensional result). Let ε, γ ∈ (0, 1]. There is an ε-differentially private algorithm A such that, for all distributions P supported in [0, 1], for all natural numbers n > polylog 1/γ ε , there exists a distribution Q (with D ∞ (P, Q) ≤ ln 2) such that the following is satisfied.For any ε-DP algorithm A ′ , with probability at least 0.75 over the randomness of x ∼ P n and additional randomness of the algorithm,W(A( Pn ), P ) ≤ polylog n • sup P ′ ∈{P,Q} R A ′ ,n ′ (P ′ ) + γ where n ′ ≈ n polylog n/γIn this one-dimensional case, our algorithm is based on DP quantile estimation. The additive γ 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 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.3 (Informal two-dimensional result). There is a polynomial time ε-differentially private algorithm A that for any distribution P on [0, 1] 2 , any integer n, and any ε-DP algorithm A ′ with probability at least 0.75, satisfies W 2 (A( Pn ), P ) ≤ (log n) O(1) supP ′ :D∞(P,P ′ )≤ln 2 E[W 2 (A ′ ( P ′ n ′ ), P ′ )],where n ′ ≈ n polylog n . Here, the expectation is taken over the internal coin tosses of A as well as over the choice of the i.i.d. samples Pn .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.4 (Informal finite metric result). Let (X , d) be an arbitrary metric space with diameter 1. There is a polynomial time ε-differentially private algorithm A such that for any distribution P on X any integer n and any ε-DP algorithm A ′ with probability at least 0.75, satisfies 1) supW(A( Pn ), P ) ≤ (log |X | • log n) O(P ′ :D∞(P,P ′ )≤ln 2 E[W(A ′ ( P ′ n ′ ), P ′ )],where n ′ ≈ n polylog n . Here, the expectation is taken over the internal coin tosses of A as well as over the choice of the i.i.d. samples Pn .Our lower bound result is actually slightly stronger than stated in Theorem 1.4 since it holds not only for ε-DP, but also for (ε, δ)-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 √ 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/ √ 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. ℓ 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/(η/ √ d)) d at an additive cost of η in the Wasserstein distance in order to create a finite metric space to apply our result on. Since the dependence on |X | in the result above is logarithmic, this translates to a d log d η multiplicative overhead term replacing the log |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 nonprivate 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 , . . . , p k }, our results lead to a rate (up to polylogarithmic factors) of i min p i (1-p i ), pi(1-pi) n + i min p i , (1 -p i ), 1εn . 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 Θ(k/εn) (for sufficiently small ε). Now, consider the following power-law distribution over support size k: p(i) ∝ 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 Õ min{ k εn , 1 √ εn } , 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.We start by describing the rate we obtain for distributions P over R.In order to state the rate, we will use q α to represent the α-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 A,n (P ) = Õ E W P, Pn + 1 εn q 1-1 εn -q 1 εn + W(P, P | q 1 εn ,q 1-1 εn ) ,The first term is E[W(P, Pn )], 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 1 εn (q 1-1 εn -q 1 ε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/εn and 1 above q 1-1/ε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 εn equally spaced quantiles of the empirical distribution differentially privately (using a known private CDF estimation algorithm), and placing roughly 1/ε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 E[W(P, Pn )] 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/εn (due to the discretization), and so we will pay 1/ε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 E.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 1 2 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 1 2 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 ∞ (P, Q) ≤ 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/εnquantile of P (scaling the density function by half) and moving it to the last 1/ε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/ε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 ε-DP algorithm to pinpoint the location of an 1 ε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 ≤ 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 + 2 i /n)f P and balance out the extra mass by setting f Q to be (1 -2 i /n)f P between [q 1-1/2 i , q 1-1/2 i+1 ). For i ≥ 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 + 1 2 )f P and balance out the extra mass by setting f Q to be (1 -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 1 √ 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 Pn by a similar term. It is more intuitive to reason about this using an alternative (equivalent) for-mulation 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 Pn in the interval [q 1/2 i+1 , q 1/2 i ) would be roughly 1 √ 2 i n (by properties of a Binomial) and hence the contribution of this interval to the area would be roughly 1 √ 2 i n q 1/2 i -q 1/2 i+1 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 Pn 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 W(P, Q) ≥ E[W(P, Pn ]. 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+ 1 √ 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 E.1.Since the main technical challenge of proving Theorem 1.4 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.5 (Hierarchically Separated Tree). A hierarchically separated tree (HST) is a rooted weighted tree such that the edges between level ℓ and ℓ -1 all have the same weight (denoted r ℓ ) and the weights are geometrically decreasing so r ℓ+1 = (1/2)r ℓ . Let D T 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 G P on the nodes of the tree where the value of a node ν is given by the weight in P of the leaf nodes in the subtree rooted at ν. For every level ℓ ∈ [D T ] of the tree, let P ℓ be the distribution induced on the nodes at level ℓ where the probability of node ν is G P (ν). Thus P ℓ is a discrete distribution on a domain of size N ℓ , where N ℓ is the number of nodes in level ℓ of the tree. Lemma 1.6 (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:W(P, Q) = 1 2 ν r ν |G P (ν) -G Q (ν)| = ℓ r ℓ TV(P ℓ , Q ℓ ),where r ν is the weight of the edge connecting ν to its parent, and the sum is over all nodes in the tree.We will call a node ν α-active under the distribution P if G P (ν) ≥ α. Let γ P (α) be the set of α-active nodes under P and γ P ℓ (α) be the set of α-active nodes at level ℓ. Then there exists an algorithm A such that given a distribution P , ε > 0, and n ∈ N,RA,n(P ) = Õ   max ℓ r ℓ x∈[N ℓ ] min P ℓ (x)(1 -P ℓ (x)), P ℓ (x)(1 -P ℓ (x)) n + x / ∈γ P ℓ (2κ) P ℓ (x) + (|γP ℓ (2κ) | -1)κ   ,where the max is over all the levels of the tree and κ = Θ( log(n) εn ). Further, this bound matches (up to logarithmic factors) the lower bound min ε-DP A ′ sup P ′ :D∞(P,P ′ )≤ln 2 E[W(A ′ ( P ′ n ′ ), P ′ )] where n ′ ≈ n polylog n . The error rate R 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 ≈ 1/(ε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 κ-active or have weight 0, so the first two terms in R N ,n,ε (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 (≈ εn) of points. While we could try to mimic the one-dimensional solution by privately estimating a solution to the ε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 κ-active nodes, then we estimate the weight of these active nodes. Since these nodes have weight greater than log(n) ε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.6. 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 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 , • • • , u k ) ∈ {±1} k , for each coordinate j ∈ [k], if u j = +1, we move Õ(κ) mass from one node in the jth pair to the other node. If u j = -1 then we leave the jth pair of nodes alone. Since each active node has mass > κ, it's clear that each resulting distribution belongs in N (P ). We can also show that these distributions form a sufficiently hard multiple hypothesis testing problem. By DP Assouad's (Lemma D.8), this allows us to lower bound the estimation error by Ω(kκ), which is within Ω of the third term when the number of active nodes is ≥ 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 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 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. 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 X , let ∆(X ) be the set of distributions on the space X . Given a logical statement a, let χ a = 0 if a is false and 1 if a is true. For example, χ 0=0 = 1 and χ 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 A.1 (D ∞ -divergence). Given two distributions P and Q with the same support, the ∞-Rényi divergence D ∞ (P, Q) = ln sup t max P (t) Q(t) , Q(t) P (t) , if P and Q are discrete, andD ∞ (P, Q) = ln sup t max f P (t) f Q (t) , f Q (t)f P (t) , if P and Q are continuous distributions on R, and have density functions. If P and Q don't have the same support, then D ∞ (P, Q) = ∞.We will use KL(P, Q) to denote the KL-divergence, H2 (P, Q) to denote the squared Hellinger divergence and TV(P, Q) to denote the total variation distance, defined later.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 A.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 asW(P, Q) = inf π E E D(t, t 0 )π(x, x 0 ),where the infimum is over all measures π on the product space E × E with marginals P and Q respectively.Finally, for one dimensional real spaces where the metric of interest is ℓ 1 norm, we will use the following equivalent formulation of Wasserstein distance extensively. Lemma A.3 (Wasserstein formula over R). Let P, Q represent probability distributions on R with finite expectation. Then, the 1-Wasserstein distance between P, Q is equal toW(P, Q) = ∞ ∞ |F P (t) -F Q (t)|dt,where the F (•) represents the cumulative distribution function.Given an metric space X , the Wasserstein metric is a well-defined metric on the set of the probability distributions over X .A number of other distances between distributions are used in this work. Definition A.4 (KL-divergence). Given two distributions P and Q with supp(P ) ⊆ supp(Q), the KL divergence KL(P, Q) = t∈supp(P ) P (t) ln P (t) Q(t) , if P and Q are discrete, andKL(P, Q) = t∈R:f P (t)>0 f P (t) ln f P (t) f Q (t)dt if P and Q are distributions on R, and have density functions. If supp(P ) ̸ ⊆ supp(Q), then KL(P, Q) = ∞. Definition A.5 (Hellinger distance). Given two distributions P and Q, the Hellinger distanceH(P, Q) = 1 √ 2 ∥ √ P - √Q∥ 2 (where we think of P and Q as vectors representing the probability masses, and the square root being component-wise.), if P and Q are discrete. If P and Q are distributions on R, and have density functions, then H(P, Q) = Definition A.6 (Total Variation distance). Given two discrete distributions P and Q, the Total Variation distance T V (P, Q) = 1 2 ∥P -Q∥ 1 , (where we think of P and Q as vectors representing the probability masses). More generally, for any two probability measures P and Q defined on (Ω, F), the total variation distance is defined as sup A∈F |P (A) -Q(A)| where P (A) represents the probability of A under measure P and likewise for Q.We use the following relationship between Hellinger distance and KL divergence. Lemma A.7. For all distributions P, Q such that KL-divergence of P, Q is well defined, we have thatH 2 (P, Q) ≤ KL(P, Q), H 2 (P, Q) ≤ T V (P, Q) (2) A.2 Differential Privacy Lemma A.8 (Post-Processing [DMNS17]). If Algorithm A : X n → Y is (ε, δ)-differentially pri- vate, and B : Y → Z is any randomized function, then the algorithm B • A is (ε, δ)-differentially private.Secondly, differential privacy is robust to adaptive composition.Lemma A.9 (Composition of (ε, δ)-differential privacy [DMNS17]). If A is an adaptive composi- tion of m differentially private algorithms A 1 , . . . , A m , where A j is (ε j , δ j ) differentially private, then A is j ε j , j δ j -differentially private.Finally, we discuss the Laplace mechanism, which we will use in one of our algorithms.Definition A.10 (ℓ 1 -Sensitivity). The ℓ 1 -sensitivity of a function f : X n → R d is ∆ f = max x,x ′ ∈X n d ham (x,x ′ )≤1 ∥f (x) -f (x ′ )∥ 1 . Lemma A.11 (Laplace Mechanism). Let f : X n → R d be a function with ℓ 1 -sensitivity ∆ f . Then the Laplacian mechanism is algorithm A f (x) = f (x) + (Z 1 , . . . , Z d ),whereZ i ∼ Lap ∆ f ε (and Z 1 , . . . , Z d are mutually independent). Algorithm A f is ε-DP.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 P ⊂ ∆(X ) on a space X and our goal is to learn the parameter θ : P → M where M is a metric space with metric d. Given an estimation algorithm A : X n → M, we can define the estimation rate3 of A to be the function R A,n : P → R + whereR A,n (P ) = E D∼P n [d(θ(P ), 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 A simply outputs the empirical mean. Then R A,n (Ber(p)) = min{p(1 -p), 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. [CL15] proposed three desiderata that a target estimation rate R n : ∆(X ) → R + should satisfy in order to be a meaningful benchmark; 1. R n (P ) varies significantly across P 2. R n is an achievable estimation rate; there exists an algorithm A and constant α such that R A,n (P ) ≤ αR n (P ) for all P ∈ P 3. Outperforming the benchmark 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 B.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.We will start by defining a target estimation rate. We'll say an algorithm is α-instance optimal if it uniformly achieves this target estimation rate up to a multiplicative α factor. For each distribution P ∈ P, we define a neighbourhood N (P ).Definition B.1. Given a function N : P → P(P), where P(P) is the power set of P, we define the optimal estimation rate with respect to N to be:EQUATIONAn algorithm A is α-instance optimal with respect to N if for all P ∈ P, R A,n (P ) ≤ αR N ,n (P )If an algorithm A uniformly achieves the optimal estimation rate wrt a function N , then this implies that for all distributions P , the error of the algorithm A on P is competitive with an algorithm that is told the additional information that the distribution is in N (P ). Given a function N , it is possible that there does not exist an algorithm that uniformly achieves R N ,n . For example, as discussed in the introduction, if N (P ) = {P }, then R N ,n is not uniformly achievable. Conversely, if N (P ) is not chosen carefully, then the estimation rate R 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 D over distributions that they care about, and their objective may be to minimizeE P ∼D [R A,n (P )].Suppose that for the A ⋆ chosen by the algorithm designer, and for our neighborhood map N , the function R A ⋆ ,n (P ) does not vary too much over N (P ) on average. Formally, let disc N A ⋆ (P ) = supP ′ ∈N (P ) (R A ⋆ ,n (P ′ ) -R A ⋆ ,n (P )and let discN A ⋆ (P ) = E P ∼D [disc N A ⋆ (P )].Then for any algorithm A that is α-instance optimal with respect to N , we can writeE P ∼D [R A,n (P )] ≤ α • E P ∼D [ sup P ′ ∈N (P ) R A ⋆ ,n (P ′ )] = α • E P ∼D [R A ⋆ ,n (P ) + disc N A ⋆ (P ))] = α • E P ∼D [R A ⋆ ,n (P )] + α • E P ∼D [disc N A ⋆ (P ))] = α • E P ∼D [R A ⋆ ,n (P )] + disc N A ⋆ (P ) .In other words, as long as the algorithm A ⋆ 's performance is relatively constant over N (P ) on average over the distribution of interest, the instance optimal algorithm (that is not tailored to D) is competitive with A ⋆ . 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 N (P ) be a small D ∞ 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 ∞ 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 Appendix D.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 instanceoptimality 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 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 instanceoptimality results with any o(K) competitive ratio. Indeed, for any two given distributions on [K] with TV distance α, Õ( 1 α 2 ) samples suffice to distinguish them, whereas learning a near uniform distribution on K atoms requires Ω(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 R, instance optimality with respect to neighborhoods of size 2 is achievable. In the standard version of this benchmark metric, N (P ) = {P, Q P } where Q P can be any distribution and is chosen to maximise R 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 ∞ ball around P (D ∞ (P, Q) ≤ ln 2), and ask that an algorithm is competitive with an algorithm that is told the additional information that P ∈ {P, Q} (in the worst case over distributions Q that are in this D ∞ ball). That is, we define the benchmark estimation rate to be R loc,n (P ) = supEQUATIONNote 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 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 √p(1-p) √ n+ min{p, 1 -p, 1 εn } which also appear to be better than the instance-optimal lower bounds in [MSU22] , which take the form√ p(1-p) √ n + 1ε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., N (P ) is {P, Q} for a worst-case Bernoulli Q. On the other hand, the notion of instanceoptimality we use would only consider Bernoullis Q such that D ∞ (P, Q) ≤ ln 2. When p is close to 0, the lower bound in [MSU22] would on this instance consider Q to be Bern(p + 1 εn ), which can have a large D ∞ -distance from P , and so isn't in the neighborhood used in our notion of instanceoptimality. 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 instanceoptimality 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 N (P ), they apply for the average case over a specific distribution over 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 R loc,n,ε to denote the version of Eqn 4 where the minimum is taken over all ε-DP mechanisms, and R N ,n,ε to define the optimal ε-DP estimation rate, i.e. Eqn 3 where the minimum is taken over all ε-DP mechanisms.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 B.2. Let M be a class of algorithms. We say that an algorithm A is α-locally minimal with respect to a neighborhood map N , if for all instance P , and allA ′ ∈ M, there is a Q ∈ N (P ) such that R A,n (Q) ≤ α • R A ′ ,n (Q).In words, local minimality says that for any other A ′ , the algorithm A is competitive with A ′ for some instance in the neighborhood of P . Put differently, no A ′ can be uniformly much better than 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. R A,n (P ) = 1; R A,n (Q) = 0; R A ⋆ ,n (P ) = 0; R A ⋆ ,n (Q) = 0;Then one can verify that A is (1-)locally minimal in M. However, it is not α-instance optimal for any finite α as it fails to satisfy the definition at P . Example B.4 (Instance Optimality ̸ ⇒ Local Minimality). Consider a set of instances {P 1 , P 2 , P 3 } with N (P 1 ) = {P 1 , P 2 }, N (P 2 ) = {P 1 , P 2 , P 3 }, N (P 3 ) = {P 2 , P 3 }. Let M contain algorithmsA, A ⋆ with R A ⋆ ,n (P 1 ) = 1; R A,n (P 1 ) = 2α; R A ⋆ ,n (P 2 ) = 2α; R A,n (P 2 ) = 4α 2 ; R A ⋆ ,n (P 3 ) = 4α 2 ; R A,n (P 3 ) = 4α 2 .Then one can verify that A is (1-)instance optimal in M. However, it is not α-locally minimal at P 1 .Under smoothness assumptions on A with respect to N , one can argue that the two notions are essentially equivalent.Proposition B.5. Let A be such that for all instances P and for allQ ∈ N (P ), R A,n (Q) ≤ β • R A,n (P ).Further, suppose that N (P ) is compact for any P . If A is α-instance optimal in M with respect to N , then it is αβ-locally minimal.Proof. Let P be an instance and let A ′ be a competing algorithm. By definition of α-instance optimality,R A,n (P ) ≤ α • sup Q∈N (P ) R A ′ ,n (Q).By compactness, this implies that there is a Q achieving the supremum. In other words, there existsQ ⋆ ∈ N (P ) such that R A,n (P ) ≤ α • R A ′ ,n (Q ⋆ ).Since Q ⋆ ∈ N (P ), our smoothness assumption implies thatR A,n (Q ⋆ ) ≤ β • R A,n (P ).Combining the last two inequalities, thisQ ⋆ satisfies R A,n (Q ⋆ ) ≤ αβ • R A ′ ,n (Q ⋆ ).Since P and A ′ were arbitrary, this implies that A is αβ-locally minimal.Proposition B.6. Let A be such that for all instances P and for allQ ∈ N (P ), R A,n (Q) ≥ β -1 • R A,n (P ). If A is α-locally minimal in M with respect to N , then it is αβ-instance optimal.Proof. Let P be an instance and let A ′ be a competing algorithm. By definition of α-local minimality, there is aQ ⋆ ∈ N (P ) such that R A,n (Q ⋆ ) ≤ α • R A ′ ,n (Q ⋆ ).Since Q ⋆ ∈ N (P ), our smoothness assumption implies thatR A,n (P ) ≤ β • R A,n (Q ⋆ ).Combining the last two inequalities, thisQ ⋆ satisfies R A,n (P ) ≤ αβ • R A ′ ,n (Q ⋆ ) ≤ αβ • sup Q∈N (P ) cost(A ′ (Q), Q).Since P and A ′ were arbitrary, this implies that A is αβ-instance optimal.A similar pair of results hold when the comparator algorithm A ′ is smooth with respect to the neighborhood map.We finish by noting relaxations of the above definitions that share the same semantic meaning (our algorithms will achieve these relaxed notions). Definition B.7. Given a function N : P → P(P), where P(P) is the power set of P, we define the optimal estimation rate with respect to N to be:EQUATIONAn algorithm A is (α, β, γ)-instance optimal with respect to N if for all P ∈ P, R A,n (P ) ≤ αR N ,βn (P ) + γ Definition B.8. Let M be a class of algorithms. We say that an algorithm A is (α, β, γ)-locally minimal with respect to a neighborhood map N , if for all instance P , and all A ′ ∈ M, there is aQ ∈ N (P ) such that R A,n (Q) ≤ α • R A ′ ,βn (Q) + γ.Note that we think of β ∈ (0, 1] and γ 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 B.1) can be obtained by setting β = 1 and γ = 0.In our work, for most settings of interest, we roughly achieve β = 1/(log n) O(1) and γ to be an arbitrarily small polynomial in the inverse of the number of samples 1/n at a log(1/γ) cost to the multiplicative factor. We don't view this as a significant issue since we expect the benchmark rate with Õ(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.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 B.1. While our results for density estimation in R 1 satisfy this notion, they also satisfy a stronger notion described in Section B.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 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 ∞ 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 Appendix B 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. [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 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 nonparametric 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.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 D.2 we'll characterise the target estimation rate R N ,n . In Section D.3, we'll then provide an ε-DP algorithm and prove that it achieves this target estimation rate up to logarithmic factors.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 D.1 (Hierarchically Separated Tree). A hierarchically separated tree (HST) is a rooted weighted tree such that the edges between level ℓ and ℓ -1 all have the same weight (denoted r ℓ ) and the weights are geometrically decreasing so r ℓ+1 = (1/2)r ℓ . Let D T 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 D.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 ∈ V• d(u, v) ≤ d T (u, v) • E[d T (u, v)] ≤ O(log M ) • d(u, v)The depth of the HST is logarithmic in the size of the metric space, D T = log M .An immediate consequence of the O(log M ) metric distortion in Lemma D.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 D.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 D.2 to translate the problem of density estimation on a bounded ball in R d into density estimation on an HST, one discretizes the metric, paying a small additive term. Corollary D.3. Given α > 0, there is a probabilistic embedding f of [0, 1] d into an HST such that for all x, y ∈ [0, 1] d :• d(x, y) -α ≤ d T (f (x), f (y)) • E[d T (f (x), f (y)] ≤ O(d • log 1 α ]) • (d(x, y) + α)The distortion is logarithmic in 1 α , so taking α 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 1 α , up to an α additive error. A distribution P on the the underlying metric space in an HST induces a function G P on the nodes of the tree where the value of a node ν is given by the weight in P of the leaf nodes in the subtree rooted at ν. For every level ℓ ∈ [D T ] of the tree, let P ℓ be the distribution induced on the nodes at level ℓ where the probability of node ν is G P (ν). Thus P ℓ is a discrete distribution on a domain of size N ℓ , where N ℓ is the number of nodes in level ℓ of the tree. Lemma D.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:W(P, Q) = 1 2 ν r ν |G P (ν) -G Q (ν)| = ℓ r ℓ TV(P ℓ , Q ℓ ),where r ν is the length of the edge connecting ν to its parent, and the sum is over all nodes in the tree.Recall the definition of our neighbourhood.N (P ) = {Q ∈ P | D ∞ (P, Q) ≤ ln 2}We will call a node ν, α-active node under the distribution P if the weight in P of the sub-tree rooted at ν is greater than α. Let γ P (α) be the set of α-active nodes under P and γ P ℓ (α) be the α-active nodes at level ℓ. Theorem D.5. Given a distribution P on [N ], ε > 0, δ ∈ [0, 1], and n ∈ N, let κ = 1 10εn min{W 0.45ε δ , 0.6} where W (x) is the Lambert W function so W (x)e W (x) = x, thenRN ,n,ε(P ) = Ω   max ℓ r ℓ x∈[N ℓ ] min P ℓ (x)(1 -P ℓ (x)), P ℓ (x)(1 -P ℓ (x)) n + x / ∈γ P ℓ (2κ) P ℓ (x) + (|γP ℓ (2κ) | -1)κ   ,where the max is over all the levels of the tree.Note that κ ≈ 1 εn min{log(1/δ), 1} so the dependence on ε and n in Theorem D.5 matches the upper bound in Theorem D.13. The error rate R N ,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 estimate the value of nodes that have probability less than 1/(ε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 1 εn -active or have weight 0, so the first two terms in R N ,n,ε (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.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 ℓ ; let P be an estimate of the distribution P and Pℓ be the induced estimate of the distribution at level ℓ. Then for any distribution PW(P, P ) = ℓ∈[D T ] r ℓ T V (P ℓ , Pℓ ).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 ℓ. Theorem D.6. For every level ℓ ∈ [D T ], define the neighborhood of P ℓ asN ℓ : ∆([N ℓ ]) → P(∆([N ℓ ])) by N ℓ (P ℓ ) = {Q ℓ | D ∞ (P ℓ , Q ℓ ) ≤ ln 2}. Then, R N ,n,ε (P ) ≥ max ℓ∈[D T ] r ℓ • R N ℓ ,n,ε (P ℓ ),where the error of P is measured in the Wasserstein distance and P ℓ is measured in the TV distance.Recall that R N ℓ ,n,ε (P ℓ ) is the optimal estimation rate with respect to N ℓ where the error is measured with respect to the total variation error. The proof of Theorem D.6 can be found in Appendix G.In light of Theorem D.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 ∈ N, and let P be a distribution on[N ]. Define N (P ) = {Q | D ∞ (P, Q) ≤ ln 2}.Our goal is to give a lower bound for R f,n,ε (P ), where the metric is the TV distance.Theorem D.7. Given ε > 0 and δ ∈ [0, 1], let κ = 1 10εn min{W 0.45ε δ , 0.6} where W (x) is the Lambert W function so W (x)e W (x) = x. Given a distribution P , R N ,n,ε (P ) = Ω   x∈[N ] min P (x)(1 -P (x)), P (x)(1 -P (x)) n +x / ∈γ P (2κ)P (x) + (|γ P (2κ) | -1)κ  Theorem D.5 follows immediately from Theorem D.6 and Theorem D.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 G for completeness.Lemma D.8. [A extension of (ε, δ)-DP Assouad's method [ASZ21] ] Let k 0 , k 1 , • • • be a sequence of natural numbers such that s k s < ∞, ε > 0 and δ ∈ [0, 1]. Given a family of distributions P ⊂ ∆(X ) on a space X , a parameter θ : P → M where M is a metric space with metric d, suppose that there exists a set V ⊂ P of distributions indexed by the product of hypercubesE k0 × E k1 × • • • where E k := {±1} k such that for a sequence τ 0 , τ 1 , • • • , ∀(u 0 , u 1 , • • • ), (v 0 , v 1 , • • • ) ∈ E k0 × E k1 × • • • , d(θ(p u ), θ(p v )) ≥ 2 s τ s ks j=1 χ u s j ̸ =v s j . (6)For each coordinate s ∈ N, j ∈ [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) = 2 |E k0 × E k1 × • • • | u∈E k 0 ×E k 1 ×•••:u s j =+1 p u , p -(s,j) = 2 |E k0 × E k1 × • • • | u∈E k 0 ×E k 1 ×•••:u s j =-1 p u ,and let ϕ s,j : X n → {-1, +1} be a binary classifier. Thenmin A is (ε,δ)-DP max p∈V R A,n (p) ≥ 1 2 s τ s ks j=1 min ϕs,j is (ε,δ)-DP ( Pr X∼p n +(s,j) (ϕ s,j (X) ̸ = 1)+ Pr X∼p n -(s,j) (ϕ s,j (X) ̸ = -1)),where the min on the LHS is over all (ε, δ)-DP mechanisms, and on the right hand side is over all (ε, δ)-DP binary classifiers. Moreover, if for all s ∈ N, j ∈ [k s ], there exists a coupling (X, Y ) between p n +(s,j) andp n -(s,j) with E[d Ham (X, Y )] ≤ D s , then min A is (ε,δ)-DP max p∈V R A,n (p) ≥ s k s τ s 2 (0.9e -10εDs -10D s δ)Note that an upper bound on T V (P i , P j ) ≤ γ implies there exists a coupling (X, Y ) between P n i and P n j such that E[d Ham (X, Y )] ≤ nγ. We will separately prove that each of the three terms in Theorem D.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 N (P ) that satisfy the conditions of Lemma D.8. Lemma D.9 and Lemma D.10 give lower bounds on the noise due to privacy. Lemma D.11 gives lower bounds based on the error due to sampling.κ = 1 10ε min{W 0.45ε δ , 0.6},where W (x) ≈ 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 ≤ κ, thenEQUATIONLemma D.9. Given a distribution P , ε > 0, δ ∈ [0, 1] and n ∈ N,R N ,n,ε (P ) ≥ 0.1 γ P 2κ n -1 κ n .Proof. Let L = |γ P (2κ/n) | 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 = ⌊L/2⌋ ≥ 1. Using the notation from Lemma D.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 ), • • • , (a + k , a - k ). Given u ∈ E k , define the distribution p u as follows: for all a b j ∈ γ P (2κ/n), P u (a + j ) = P (a + j )+(κ/n) and P u (a - j ) = P (a - j )-(κ/n) if u j = +1 and P u (a + j ) = P (a + j )-(κ/n) and P u (a - j ) = P (a - j )+(κ/n) if u j = -1. For all other x, P u (x) = P (x). It is immediate that for all u, P u ∈ N (P ). For any pair u, v, d(θ(p u ), θ(p v )) = TV(p u , p v ) = d Ham (u, v)(κ/n), so that Equation ( 6) is satisfied with τ = 1 2 (κ/n). Further, given j ∈ [k], p +j and p -j only differ on the probability of a + j and a - j , so D/n = max j TV(p +j , p -j ) = κ/n and by Equation (7), 0.9e -10εD -10Dδ ≥ 0.2. Noting that k ≥ (1/2)(γ P (2κ/n) -1) completes the proof.Lemma D.10. For all ε > 0, δ ∈ [0, 1], n ∈ N and distributions P on[N ], if κ < n/2, then R N ,n,ε (P ) ≥ Ω   x / ∈γ P (2κ) P (x)   . Since κ ≤ 1εn , the condition that κ < n/2 is a mild condition. For example, it is satisfied whenever ε > 2/n 2 . Similar to the proof of Lemma D.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 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 ∈ N, letS s = {x ∈ [N ] | P (x) ∈ (2 -s-1 , 2 -s ]}. Proof. Given s ∈ N, let S ′ s = S s ∩ {x | P (x) ≤ 2κ/n} and d s = |S ′ s |.Let us first consider the case that there exists a scale s * with d s * = 1 and P (x * ) ≥ 1 8x / ∈γ P (2κ) P (x) where x * is the element in S ′ s * . Define P ′ by P ′ (x * ) = (1/2)P (x * ) and for all x ̸ = x * , P ′ (x) = 1-(1/2)P (x * ) 1-P (x * ) P (x). Since (1/2)P (x * ) ≤ 2κ/n ≤ 1/2, P ′ ∈ N (P ). In this case we will use Lemma D.8 with k 0 = 1, k s = 0 otherwise, and E k0 corresponds to the set of distributions {P, P ′ }. Then noting that TV(P, P ′ ) = (1/2)P (x * ) and using eqn (7) we have that τ = 1 4 P (x * ) and D = (1/2)P (x * )n ≤ κ so that R N ,n,ε (P ) ≥ (1/8)P (x * )(0.2) = Ωx / ∈γ P (2κ) P (x) and we are done.Next suppose that for all scales s such that d s = 1 we have P (x * ) < 1 8x / ∈γ P (2κ) P (x). Let s * be the smallest s such that d s = 1. Since the scales 2 -s-1 are geometrically decreasing,s:ds=1 x∈Ss∩{x|P (x)≤2κ/n} P (x) ≤ 2 s:ds=1 2 -s-1 ≤ 4 • 2 -s * -1 ≤ 1 2x / ∈γ P (2κ)It follows that s:ds>1 x∈Ss∩{x|P (x)<2κ/n} P (x) ≥ (1/2) x / ∈γ P (2κ) P (x). Further,x / ∈γ P (2κ)P (x) ≤ 2 s:ds>1 x∈Ss∩{x|P (x)<2κ/n} P (x) ≤ 4 s:ds>1 2 -s-1 d s ≤ 16 s:ds>1 2 -s-1 ⌊d s /2⌋.Thus we can (up to constants) ignore scales such that d s ≤ 1 and assume that d s is even for all scales.Let k s = ⌊d s /2⌋. Now, within each scale S s , pair the elements to form k s distinct pairs (a + s,j , a - s,j ).Given (u 0 , u 1 , • • • ) ∈ E k0 × E k1 × • • • , 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 ∈ N (P ). Further, using notation from Lemma D.8, Equation (6) is satisfied with τ s = 1 2 2 -s-2 since d(θ(p u ), θ(p v )) = TV(p u , p v ) = 2 s τ s d Ham (u s , v s ) and D s /n = max j TV(p +(s,j) , p -(s,j) ) = 2 -s-2 , which is less than κ whenever k s > 0. By eqn (7), 0.9e -10εDs -10D s δ ≥ 0.2 whenever k s > 0 so by Lemma D.8 we haveR N ,n,ε (P ) ≥ s 1 2 2 -s-2 k s (0.2) = 0.2 4 s:ds>1 2 -s-1 ⌊d s /2⌋ ≥ 0.2 4 × 16 x / ∈γ P (2κ) P (x)which completes the proof.Next we lower bound the statistical term. Lemma D.11. For all n ∈ N, ε > 0, δ ∈ [0, 1] and distributions P , if n ≥ 2 and ε > 2/n, thenR N ,n,ε (P ) ≥ R N ,n (P ) ≥ Ω   x∈[N ] min P (x)(1 -P (x)), P (x)(1 -P (x) n   .To streamline the notation, we will use L(x) to denote min x(1 -x), x(1-x) n. In order to prove Lemma D.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 D.11. Lemma D.12. Given any pair of distributions P and Q on the same domain,min ϕ Pr X∼P n (ϕ(X) = 1) + Pr X∼Q n (ϕ(X) = -1) ≥ 1 2 (1 -nKL(P, Q)),where the minimum is over all binary classifiers. In particular, if P = Bernoulli(p -α) andQ = Bernoulli(p + α) where 0 ≤ α ≤ 1 2 L(p) then min ϕ Pr X∼P n (ϕ(X) = 1) + Pr X∼Q n (ϕ(X) = -1) ≥ 1/4,where again the minimum is over all binary classifiers.The proof of Lemma D.12 can be found in Appendix G Proof of Lemma D.11. As in the proof of Lemma D.10, first suppose there exists a scale s * with d s * = 1 and there exists x * ∈ S s * such that1 2 L(P (x * )) ≥ 1 60 x∈[N ]L(P (x)).Then define a distribution P ′ by P ′ (x * ) = P (x * ) -1 2 L(P (x * )) and for all x ̸ = x * , P ′ (x) = 1-P (x * )+ 1 2 L(P (x * )) 1-P (x * ) P (x). Then P ′ ∈ N (P ) since 1 2 L(P (x * )) < 1 2 min{P (x * ), (1 -P (x * ))}. Then we will use Lemma D.8 with k 0 = 1 and k s = 0 for s > 0, and E k0 corresponds to {P, P ′ }. Now, KL(P ′ , P ) = (P (x * ) -1 2 L(P (x * ))) ln P (x * ) -1 2 L(P (x * )) P (x * ) + (1 -P (x * ) + 1 2 L(P (x * ))) ln 1 -P (x * ) + 1 2 L(P (x * )) 1 -P (x * ) ≤ 1 4n(for more detail on the proof of this inequality see the proof of Lemma D.12) somin ϕ Pr X∼P n (ϕ(X) = 1) + Pr X∼P ′ n (ϕ(X) = -1) ≥ 1 2 (1 -nKL(P, P ′ )) ≥ 1/4and τ 0 = TV(P, P ′ ) = 1 2 L(P (x * )). Thus by Lemma D.8,R N ,n,ε (P ) ≥ R N ,n (P ) ≥ 1 2 L(P (x * )) 1 4 ≥ 1 480 x∈[N ] 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 haveL(P (x s )) ≤ 1 30 x∈[N ] L(P (x)),where S s = {x s }. As in the proof of Lemma D.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 s:ds=1L(P (x s )) ≤ χ d0=1 L(P (x 0 )) + s>0:ds=1 min 2 -s , 2 -s n ≤ χ d0=1 L(P (x 0 )) + (2 + √ 2) min 2 -s * , 2 -s * n ≤ χ d0=1 L(P (x 0 )) + 2(2 + √ 2) min P (x s * ), P (x s * ) n ≤ χ d0=1 L(P (x 0 )) + 4(2 + √ 2)L(P (x s * )) ≤ 1 + 4(2 + √ 2) 30 x∈[N ] L(P (x)).Therefore, s:ds>1 x∈Ss L(P (x)) ≥ (1/2) x∈[N ] L(P (x)) and soEQUATIONwhere the first inequality follows from ⌊d s /2⌋ ≥ (1/3)d s whenever d s > 1, and the second follows because 2 -s-1 ≤ P (x) ≤ 1/2 for all x ∈ S s such that d s > 1.Assume that d s is even for all s. Within each scale 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 ∈ N, let α s = 1 2 L(2 -s-1 ), and note that for all x ∈ S s and s > 0, α s ≤ 1 2 L(P (x)). Given (u 0 , u 1 , • • • ) ∈ E k0 ×E k1 ×• • • , define p u by p u (a + s,j ) = p(a + s,j )+α s and p u (a - s,j ) = p(a - s,j ) -α s if u s j = +1 and p u (a + s,j ) = p(a + s,j ) -α s and p u (a - s,j ) = p(a - s,j ) + α s if u s j = -1. For all other elements, p u (x) = p(x). Then, for all u, p u ∈ N (P ). Further, using notation from Lemma D.8, we have τ s = α 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 ) + α s and p -(s,j) (a + s,j ) = P (a + s,j ) -α s . Therefore, by Lemma D.12, and the post-processing inequality,min ϕ Pr X∼p +(s,j n (ϕ(X) = 1) + Pr X∼p -(s,j n (ϕ(X) = -1) ≥ 1/4.Lemma D.8 then implies the result. Theorem D.7 follows immediately from Lemma D.9, Lemma D.10 and Lemma D.11.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 (≈ εn) of points. While we could try to mimic the onedimensional solution by privately estimating a solution to the ε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 log(1/δ) εn -active nodes, then we estimate the weight of these active nodes. Since these nodes have weight greater than log(1/δ) εn , we can privately estimate them to within constant multiplicative error.Let X be the underlying metric space so P ∈ ∆(X ). For any set S of nodes and a function F defined on the nodes, define the function F | S as F | S (ν) = F (ν) if ν ∈ S and F | S (ν) = 0 otherwise. Given two functions F and G defined on the nodes, we defineW(F, G) = ν r ν |F (ν) -G(ν)|,where r ν is the length of the edge connecting ν to its parent, and the sum is over all nodes in the tree. So by Lemma D.4, W(P, Q) = W(G P , G Q ). Note that W satisfies the triangle inequality.Algorithm 1 PrivDensityEstTree1: Input: D ∈ X n , ε 2: G P = EmpDist(D)▷ Compute empirical distribution. G P (ν) + |γ P ℓ 1 2εn -1| log(1/β) εnThis bound has the same three terms as our lower bound on R N ,n,ε in Theorem D.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 ≈ log(1/δ)/ε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 D T of the lower bound. Since we can not hope to locate the set of log(1/δ)/(εn)-active nodes exactly with a private algorithm, we find a set γn that is guaranteed to satisfyγ P max 2 εn + 2 log(2/β) n , 192 log(n/β) n ⊂ γε ⊂ γ P 1 2εn . Note that max 2 εn + 2 log(2/β) n , 192 log(n/β) n ≤ C log(n/β)εn so the error introduced here by not estimating γ P 1 εn perfectly is at most a logarithmic multiplicative factor. 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), P (x)/n} even when P (x) is small. The proof of the following lemma is contained in Appendix G. Lemma D.14. For any distribution P , if log(n/β) > 1 then with probability 1 -3D T β,W( G P , G P ) ≤ ℓ∈[D T ] x∈[N ℓ ] min P ℓ (x)(1 -P ℓ (x)), 4 3 P ℓ (x)(1 -P ℓ (x)) log(n/β) n Algorithm 2 EmpDist 1: Input: D ∈ X n , A 2:Let Pn be the empirical distribution.3: for all node ν do 4: for all ν ∈ γε,ℓ doG P (ν) =        0 G Pn (ν) < √ log(n/β) n 1 G Pn (ν) > 1 - √ log n/β) n G Pn (ν) otherwise6:for all children ν ′ of ν do7: if G P (ν ′ ) + Lap( 1 εn ) > 2κ + log(2/β) εn then 8:γε,ℓ+1 = γε,ℓ+1 + {ν ′ } 9:ℓ = ℓ + 1 10: return ∪γ ε,ℓThe goal of Algorithm 3 is to estimate the set of 1/(εn)-active nodes.The next lemma allows us to bound how close to the goal we get. The proof is contained in Appendix G. Lemma D.15. Let γε be the set of active nodes found in Algorithm 1. Then with probability 1 -D T (log n + 4εn)β,γ P max 2 εn + 4 log(2/β) εn , 192 log(n/β) n ⊂ γε ⊂ γ P 1 2εn.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 D.16.If γ P max{ 2 εn + 4 log(2/β) εn , 192 log(n/β) n } ⊂ γε then W( G P , G P | γε ) ≤ W(G P , G P ) + W(G P , G P | γ P (max{ 2 εn +4 log(2/β) εn , 192 log(n/β) n }) )The key component of this proof is that any discrepancy between the weight of the nodes on P and that assigned by G P was already paid for in W(P, G P ). The final step in Algorithm 1 is to project the noisy function G Pn, γε 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 W distance between P and G Pn, γε . 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 D.17. For any real-valued function G on the nodes of the HST such that G(ν 0 ) = 1 where ν 0 is the root node and given any distribution P ,W(P, Projection(G)) ≤ 4W(G P , G).Combining the above lemmas appropriately gives the proof of Theorem D.13 (see Appendix G). for all nodes ν at level ℓ do 5:Let A ν = G(ν ′ )where the sum is over the children of ν.Let d ν be the number of children of ν 7:if A ν = 0 then 8:for all children ν ′ of ν do 9:Ḡ(ν ′ ) = 1 dν Ḡ(ν) 10: else 11:for all children ν ′ of ν do 12:Ḡ(ν ′ ) = Ḡ(ν) Aν G(ν ′ ) 13: return ḠProof of Theorem D.13. The privacy follows from the fact that each user contributes to at most D T queries in LocateActiveNodes and at most one coordinate in the computation of G Pn,γn in line 4 in PrivDensityEstTree.For the utility, we will consider each level ℓ individually. First suppose that |γ P ℓ (1/(2εn)) | > 1.W(P ℓ , ( Pε ) ℓ ) ≤ 2W((G P ) ℓ , ( G Pn,γε ) ℓ ) ≤ 2 W((G P ) ℓ , ( G P ) ℓ ) + W(( G P ) ℓ , ( G P |γ ε ) ℓ ) + W(( G P |γ ε ) ℓ , ( G Pn,γε ) ℓ ) (9) ≤ 2 2W((G P ) ℓ , ( G P ) ℓ ) + W((G P ) ℓ , (G P | γ P (max{ 2 εn +2 log(2/β) n , log(n/β) n }) ) ℓ ) + W(( G P | γε ) ℓ , ( G Pn,γε ) ℓ )where the first inequality follow from Lemma D.17, the second inequality follows from the triangle inequality and Lemma D.4, and the third follows from Lemma D.16 and Lemma D.15. Finally,W( G P | γε ) ℓ , ( G Pn,γε ) ℓ ) ≤ ν∈γ P ( 1 2εn ) r ν |Lap( 1 εn )| ≤ 1 2 ν∈γ P ( 1 2εn ) r ν |Lap( 1 εn )|The final statement then follows from Lemma D.14 and basic concentration bounds on the Laplacian distribution.If |γ P ℓ (1/(2εn)) | = 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 ∈ γ P ℓ (1/(2εn)). First suppose that P ℓ (x) > 1 -1 2εn then no node that is in a level above x, but is not a direct ancestor of x is in γ P ℓ (1/(2εn)). Therefore, since the projection algorithm is top-down, ( Pn,ε ) ℓ will be concentrated on x. Therefore, the error of level ℓ 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 ℓ (x) < 1 -1 2εn then sum of the inactive nodes (in term two) dominates the error due to adding noise to P (x)Let us now consider the setting of estimating distributions P on X = 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 ∞ (P, Q P ) ≤ 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 Appendix B.1 and Appendix B.3. In Appendix E.1, we lower bound this estimation rate using hypothesis testing techniques. Then, in Appendix E.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 J.To state the main theorem in this section, we will introduce some notation. We start by defining the restriction of a distribution. Definition E.1. For any distribution P over R with a density function, the restriction P | u,v of P with respect to u ≤ v ∈ R is defined as the distribution with the following CDF function F':F ′ Pu,v (t) =    0 t < u F P (x) u ≤ t < v 1 t ≥ vIf u = v, then F ′ 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 E.2. For 0 < α ≤ 1, the α-quantile of a distribution P over R is defined as follows:q α (P ) = arg min t { Pr y∼P (y ≤ t) ≥ α}.When the distribution P is clear from context, we will sometimes abuse notation and use q α when we mean q α (P ). The main theorem we will prove in this section is the following: Theorem E.3. There exists a constant C such that given a continuous distribution P on R with bounded expectation and ε∈ (0, 1], n ∈ N, R loc,n,ε (P ) = Ω 1 εn q 1-1 Cεn -q 1 Cεn + W(P, P | q 1 Cεn ,q 1-1 Cεn ) + 1 √ log n E W(P | q 1 Cεn ,q 1-1 Cεn , Pn | q 1 Cεn ,q 1-1 Cεn ) ,where Pn is the empirical distribution on n samples drawn independently from P .The same result can be extended to (ε, δ)-DP algorithms as well for δ = o( 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 C), 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 distributionlonger 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 E.1.1 we characterize the cost of private instance optimality, and in Section E.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.The main theorem we will prove in this section is the following. Theorem E.4. Fix ε ∈ (0, 1], n ∈ N. For all distributions P over R that have a density function and finite expectation, there exists another distribution Q ′′ such that D ∞ (P, Q) ≤ 2, that is indistinguishable from P given O(n) samples such that for all ε-DP algorithms A : R n → ∆(R), with probability at least 0.25 over the draws x ∼ P n , x ′ ∼ Q ′′n , the following holds for some constant C.A(x)), W(Q ′′ , A(x ′ ))) ≥ 1 4Cεn q 1-1 Cεn -q 1 Cεn + 1 4 W(P, P | q 1 Cεn ,q 1-1).We start with some notation. For any distribution P with a density, let f P denote its density function.Throughout this section, we will use q α to represent the α-quantile of distribution P . Let L(P ) be the 'starting point' of distribution P (defined as inf t∈R {t : F P (t) > 0} if the infimum exists, and -∞ otherwise.Next, we describe some results on differentially private testing that we will use. We say that a testing algorithm A test distinguishes two distributions P and Q with n samples, if given the promise that a dataset of size n is drawn from either P n or Q n , with probability at least 2 3 , it outputs P if the dataset was drawn from P n and Q if it was drawn from Q n . We now state a theorem lower bounding the sample complexity of differentially private hypothesis testing. Theorem E.5 ([CKM + 19, Theorem 1.2]). Fix n ∈ N, ε > 0. For every pair of distributions P, Q over R, if there exists an ε-DP testing algorithm4 A test that distinguishes P and Q with n samples, thenn = Ω 1 ετ (P, Q) + (1 -τ (P, Q))H 2 (P ′ , Q ′ ),whereτ (P, Q) = max R max{e ε f P (t) -f Q (t), 0}dt, R max{e ε f Q (t) -f P (t), 0}dt ,and H 2 (•, •) is the squared Hellinger distance between P ′ = min(e ε Q,P ) 1-τ (P,Q) , andQ ′ = min(e ε ′ P,Q) 1-τ (P,Q) , where 0 ≤ ε ′ ≤ ε is such that if τ (P, Q) = R max{f P (t) -e ε f Q (t), 0}dt, then ε ′ is the maximum value such that τ (P, Q) = R max{f Q (t) -e ε ′ f P (t), 0}dt, else ε ′ is the maximum value such that τ (P, Q) = R max{f P (t) -e ε ′ f Q (t), 0}dt.We now are ready to start proving our main theorem.Proof. (of Theorem E.4) The idea is to construct Q from P by moving mass from the leftmost quantiles to the rightmost quantile. We do this such that Q is statistically close enough to P such that the two distributions can not be distinguished with n samples, but is also far from P in Wasserstein distance. This produces a lower bound of (1/2)W(P, Q) on how well an algorithm can simultaneously estimate P and Q since if there was an algorithm that produced good estimates of P and Q in Wasserstein distance with n samples, then we could tell them apart, and this would give a contradiction.Let k be a quantity to be set later. Formally, we define Q as the distribution with the following density function.f Q (t) =    1 2 f P (t), for t < q 1/k f P (t), for q 1/k ≤ t < q 1-1 k 3 2 f P (t) for q 1-1 k ≤ t   Note that by the definition of Q, we have that D ∞ (P, Q) ≤ 2.We will prove that the sample complexity of telling apart P and Q under (ε, δ)-DP is Ω(k/ε), using known results on hypothesis testing. Then, we will argue that the Wasserstein distance between P and Q is sufficiently large. Setting k appropriately will complete the proof.Define SC ε,δ (P, Q) to be the smallest n such that there exists an (ε, δ)-DP testing algorithm that distinguishes P and Q; called the sample complexity of privately distinguishing P and Q.Lemma E.6. SC ε,δ (P, Q) = Ω(k/ε).The proof of this lemma is in Appendix J. We next argue that P and Q are sufficiently far away in Wasserstein distance.Lemma E.7.W(P, Q) ≥ 1 2k (q 1-1 k -q 1/k ) + 1 2 W(P, P | q 1 k ,q 1-1 k ).The proof of this lemma is also in Appendix J.Finally, we are ready to prove the theorem. Assume that with probability larger than 0.75 over the draw of two datasets x ∼ P n , x ′ ∼ Q n , and the randomness used by invocations of algorithm A we have that max(W(P, A(x)), W(Q, A(x ′ )) < 1 2 W(P, Q). Then, given a dataset x ′′ of size n, we can perform the following test: run the differentially private algorithm A on the dataset x ′′ and compute W(P, A(x ′′ )) and W(Q, A(x ′′ )) and output the distribution with lower distance. Then, note that W(P, Q) ≤ W(P, A(x ′′ )) + W(Q, A(x ′′ )) which implies that with probability at least 0.75, W(Q, A(x ′′ )) > 1 2 W(P, Q) if the dataset x ′′ was sampled from P n (by the accuracy guarantee). A similar argument shows that with probability at least 0.75, W(P, A(x ′′ )) > 1 2 W(P, Q) if the dataset x ′′ was sampled from Q n . Hence, with n samples we have defined a test that distinguishes P and Q. However, for k = Cεn for some constant C, by Lemma E.6 we get that any differentially private test distinguishing P and Q requires more than n samples, which is a contradiction. Hence, with probability at least 0.25 over the draw of two datasets x ∼ P n , x ′ ∼ Q n , and the randomness used by invocations of algorithm A we have that max(W(P,A(x)), W(Q, A(x ′ )) ≥ 1 2 W(P, Q) ≥ 1 4Cεn (q 1-1 Cεn -q 1/Cεn ) + 1 4 W(P, P | q 1 Cεn ,q 1-1) ,where the last inequality is by invoking Lemma E.7 with k = Cεn.In this section, we prove the following result. Theorem E.8. Fix sufficiently large natural numbers n, k > 0 and let C, C ′ > 0 be sufficiently small constants. For all algorithms A : R n → ∆ R , the following holds. For all continuous distributions P over R with a density and with bounded expectation, there exists another distribution Q (with D ∞ (P, Q) ≤ ln 2), that is indistinguishable from P given O(n) samples, such that with probability at least 0.25 over the draws x ∼ P n , x ′ ∼ Q n , the following holds.A(x)), W(Q, A(x ′ ))) ≥ C ′ √ log n E x ′′ ∼P n W P | q 1 k ,q 1-1 k , Pn | q 1 k ,q 1-1 k ,where q α is the α-quantile of P .Before going into the proof, we state the following result on the sample complexity of testing. This is a folklore result but for a proof of the lower bound see [BY02] and the upper bound see [Can17] . Theorem E.9. Fix n ∈ N, ε > 0. For every pair of distributions P, Q over R, if there exists a testing algorithm A test that distinguishes P and Q with n samples, thenn = Ω 1 H 2 (P, Q) ,wherer H 2 (•, •) represents the squared Hellinger distance between P and Q.Throughout the proof, we will use q α to represent the α-quantile of distribution P .Proof of Theorem E.8. Q is constructed by adding progressively more mass to P up until q 1/2 and subtracting proportionate amounts of mass from P afterwards. Intuitively, this is done in such a way that to 'change' P to Q, for all i ≥ 2 one has to move roughly min{ 1 √ 2 i n , 1 2 i } mass from q 1/2 i to q 1-1/2 i . This ensures that the Wasserstein distance between P and Q is larger than the expected Wasserstein distance between P and its empirical distribution on n samples Pn . This is carefully done to ensure that P is indistinguishable from Q.Formally, consider i in the range [2, log n -1). For all x ∈ (q 1/2 i , q 1/2 i-1 ], we set fQ (t) = f P (t) 1 + 2 i n . For all t ∈ (q 1-1/2 i-1 , q 1-1/2 i ], we set f Q (t) = f P (t) 1 -2 i n . Next, consider i in the range [log n, ∞). For all t ∈ (q 1/2 i , q 1/2 i-1 ], we set f Q (t) = f P (t) 1 + 1 2 . For all t ∈ (q 1-1/2 i-1 , q 1-1/2 i ], we set f Q (t) = f P (t) 1 -1 2 .Note that P has bounded expectation by assumption, and hence, so does Q. Additionally, note that D ∞ (P, Q) ≤ ln 2.There are two key considerations balanced in the design of Q. On one hand, we need Q to be indistinguishable from P given Õ(n) samples. On the other hand, we need Q to be sufficiently far away from P in Wasserstein distance. This ensures that given an accurate algorithm for estimating the density of the distribution (in Wasserstein distance) given access to Õ(n) samples from it, we can design a test distinguishing P and Q with that many samples, thereby contradicting their indistinguishability.Detailed proofs of claims below can be found in Appendix J. First, we show that P is indistinguishable from Q.Lemma E.10.KL(P, Q) = O(log n/n).Next, we establish a lower bound on the Wasserstein distance between P and Q.Lemma E.11.W(P, Q) ≥ 1 4   log n-1 j=2 1 √ 2 j n q 1-1/2 j -q 1/2 j + ∞ j=log n 1 2 j q 1-1/2 j -q 1/2 j   .Next, we upper bound the expected Wasserstein distance between the distribution P and its empirical distribution on n samples.Lemma E.12.E[W(P, Pn )] ≤ 8   log n-1 i=2 1 √ 2 i n q 1-1/2 i -q 1/2 i + ∞ i=log n 1 2 i q 1-1/2 i -q 1/2 i  We now prove a simple claim regarding restrictions.Claim E.13 (Restrictions preserve Wasserstein distance). For all datasets x, and any natural number k > 1 we have thatW(P | q 1 k ,q 1-1 k , Pn | q 1 k ,q 1-1 k ) ≤ W(P, Pn ).Finally, we are ready to put the above lemmas together to prove Theorem E.8. Fix n ′ = n C log n . Assume, for sake of contradiction, that with probability larger than 0.75 over the draw of two datasetsx ∼ P n ′ , x ′ ∼ Q n ′, and the randomness used by invocations of algorithm A we have that max(W(P,A(x)), W(Q, A(x ′ )) ≤ 1 2 W 1 (P, Q).Then, given a dataset x ′′ of size n ′ , we perform the following test: run the differentially private algorithm A on the dataset x ′′ and compute W(P, A(x ′′ )) and W(Q, A(x ′′ )) and output the distribution with lower distance. Then, note that W(P, Q) ≤ W(P, A(x ′′ )) + W(Q, A(x ′′ )) which implies that with probability at least 0.75,W(Q, A(x ′′ )) ≥ 1 2 W(P, Q) if x ′′ ∼ P n ′(by the accuracy guarantee). A similar argument shows that with probability at least 0.75,W(P, A(x ′′ )) ≥ 1 2 W(P, Q) if x ′′ ∼ Q n ′. Hence, with n ′ samples we have defined a test that distinguishes P and Q. However, by Lemma E.10 bounding the KL divergence between P and Q, Theorem E.9 on sample complexity lower bounds for testing, and Lemma A.7 on the relationship between KL and Hellinger distance, we get that any statistical test distinguishing P and Q requires more than n ′ samples, which is a contradiction. Hence, with probability at least 0.25 over the draw of two datasetsx ∼ P n ′ , x ′ ∼ Q n ′, and the randomness used by invocations of algorithm A we must have thatEQUATIONNext, note that by Lemma E.12 (with value n ′ ), we have thatE[W(P, Pn ′ )] ≤ 8   log n ′ -1 i=2 1 √ 2 i n ′ q 1-1/2 i -q 1/2 i + ∞ i=log n ′ 1 2 i q 1-1/2 i -q 1/2 i   = 8   log n C log n -1 i=2 √ C log n √ 2 i n q 1-1/2 i -q 1/2 i + ∞ i=log n C log n 1 2 i q 1-1/2 i -q 1/2 i   = 8 C log n log n-log(C log n)-1 i=2 1 √ 2 i n q 1-1/2 i -q 1/2 i + log n-1 i=log n-log(C log n) 1 2 i q 1-1/2 i -q 1/2 i + ∞ i=log n 1 2 i q 1-1/2 i -q 1/2 iAnalyzing the middle term in the above sum, we have thatlog n-1 i=log n-log(C log n) 1 2 i q 1-1/2 i -q 1/2 i ≤ log n-1 i=log n-log(C log n) 1 √ 2 i 1 √ 2 log n-log(C log n) q 1-1/2 i -q 1/2 i ≤ log n-1 i=log n-log(C log n) 1 √ 2 i √ C log n √ n q 1-1/2 i -q 1/2 i = C log n log n-1 i=log n-log(C log n) 1 √ 2 i n q 1-1/2 i -q 1/2 iSubstituting this back in the previous sum, we have thatE[W(P, Pn ′ )] ≤ 8 C log n log n-log(C log n)-1 i=2 1 √ 2 i n q 1-1/2 i -q 1/2 i + C log n log n-1 i=log n-log(C log n) 1 √ 2 i n q 1-1/2 i -q 1/2 i + ∞ i=log n 1 2 i q 1-1/2 i -q 1/2 i ≤ 8 C log n log n-1 i=2 1 √ 2 i n q 1-1/2 i -q 1/2 i + ∞ i=log n 1 2 i q 1-1/2 i -q 1/2 i ≤ 16 C log n ′ log n-1 i=2 1 √ 2 i n q 1-1/2 i -q 1/2 i + ∞ i=log n 1 2 i q 1-1/2 i -q 1/2 iwhere in the last inequality we use the fact that n ′ ≥ √ n. Hence, by Lemma E.11 (which gives a lower bound on W(P, Q)) in conjunction with the above equation, we have thatW(P, Q) ≥ C ′ √ log n ′ E[P, Pn ′ ]for some sufficiently small constant C ′ . Substituting back in Equation 10, we have that with probability at least 0.25 over the draw of two datasetsx ∼ P n ′ , x ′ ∼ Q n ′, and the randomness used by invocations of algorithm A we have thatmax(W(P, A(x)), W(Q, A(x ′ )) ≥ 1 2 C ′ √ log n ′ E[W(P, Pn ′ )] ≥ 1 2 C ′ √ log n ′ E W(P | q 1 k ,q 1-1 k , Pn ′ | q 1 k ,q 1-1 k ) ,as required.In this section, we describe an algorithm that achieves the instance optimal rate described in the previous section (up to polylogarithmic factors in some of the terms).We will be looking at distributions P supported on a discrete, ordered interval {a, a+γ, . . . , b-γ, b}. Note that by a simple coupling argument, any continuous distribution P cont on [a, b] is at most γ away in Wasserstein distance from a distribution on this grid. The dependence on γ in our bounds for discrete distributions will be inverse polylogarithmic (or better), and so our algorithms for estimating distributions P in the interval {a, a+γ, . . . , b-γ, b} also work to give similar bounds for continuous distributions on [a, b], up to a small additive factor of γ, which can be set to any inverse polynomial in the dataset size without significantly affecting our bounds.Formally, we will prove the following theorem (See Theorem E.15 for a more detailed statement).Theorem E.14. Fix ε, β ∈ (0, 1], a, b ∈ R, and γ < b -a ∈ R such that b-a γ is an integer. Let n ∈ N > c 2 log 4 b-a βγ εfor some sufficiently large constant c 2 . There exists an ε-DP algorithm A that for any distribution P on {a, a + γ, a + 2γ, . . . , b -γ, b} satisfies the following. When run with input a random sample x ∼ P n , A outputs a distribution P DP such that with probability at least 1 -β over the randomness of x and the algorithm,W(P, P DP ) = O 1 k q 1-1 k -q 1 k + W(P, P | q 1 k ,q 1-1 k ) + log n β E W P | q 1 k ,q 1-1 k , Pn | q 1 k ,q 1-1 k ,where Pn is the empirical distribution on n samples drawn independently from P , q α represents the α-quantile of distribution P , and k = ⌈ εn 4c3 log 3 b-a βγ log n β ⌉ for a sufficiently large constant c 3 .Since k ≈ εn/ log(n), this upper bound matches the lower bound in Theorem E.3 in its dependence on ε and its dependence on n (up to logarithmic factors in n). The algorithm that we will analyze proceeds by estimating sufficiently many quantiles from the empirical distribution and distributing mass evenly between the chosen quantiles. The number of quantiles is chosen carefully to ensure that the estimated α-quantiles are also approximately α-quantiles for the empirical distribution (and hence also approximately for the true distribution), and to ensure that the CDF of the output distribution closely tracks the CDF of the empirical distribution. Through a careful analysis, we are able to leverage these properties to give instance optimality guarantees for the accuracy of the algorithm.Algorithm 5 is our algorithm for density estimation, and proceeds by differentially privately estimating sufficiently many quantiles of the distribution and placing equal mass on each of them. We argue that a simple CDF based differentially private quantiles estimator A quant satisfies a specific guarantee that will be key to our analysis. See Appendix I for more details about the quantiles algorithm and formal statements and proofs therein.Algorithm 5 Algorithm A for estimating a distribution on R Input: x = (x 1 , . . . , x n ) ∼ P n , privacy parameter ε, interval end-points a, b, granularity γ, access to algorithm A quant Output:Distribution P DP on R. 1: Let k be set to ⌈ εn 4c3 log 3 b-a βγ log n β⌉ for a sufficiently large constant c 3 .2: Use Algorithm A quant referenced in Theorem I.2 with inputs interval end points a, b, granularity γ, x = (x 1 , . . . , x n ) ∈ {a, a + γ, . . . , b -γ, b} n , and desired quantile values α = {1/2k, 3/2k, 5/2k, . . . , (2k -1)/2k}, and let the outputs be q1 . . . , qk . 3: for j ∈ [k] do 4:Set P DP (q j ) = 1 k . 5: Output P DP .Observe that Algorithm 5 inherits the privacy of A quant , since it simply postprocesses the quantiles it receives from that subroutine, and hence is also ε-DP. Now, we are in a position to state our main theorem, which bounds the Wasserstein distance between the distribution output by our algorithm, and the underlying probability distribution P .Theorem E.15. Fix ε, β ∈ (0, 1], a, b ∈ R, and γ < b -a ∈ R such that b-a γ is an integer. Let n ∈ N > c 2 log 4 b-a γβε εfor some sufficiently large constant c 2 . Let P be any distribution supported on {a, a + γ, a + 2γ, . . . , b -γ, b}, and x ∼ P n . Then, Algorithm 5, when given inputs x, privacy parameter ε, interval end points a, b, and granularity γ, outputs a distribution P DP such that with probability at least 1 -O(β) over the randomness of x and the algorithm,W(P, P DP ) ≤ c log n•E W(P | q 1 k ,q 1-1 k , Pn | q 1 k ,q 1-1 k +C ′′ W(P, P | q 1 k ,q 1-1 k )+ 2 k q 1-1/k -q 1/k ,where Pn is the uniform distribution on x, q α represents the α-quantile of distribution P , c, C ′′ are sufficiently large constants, and k = ⌈εn 4c3 log 3 b-a βγ log n β ⌉, where c 3 is a sufficiently large constant.We note that using more sophisticated differentially private CDF estimators to estimate quantiles (such as ones in [BNSV15, CLN + 23]), we can also obtain a version of the same theorem for approximate differential privacy, with a better dependence on the size of the domain b-a γ (only log * ( b-a γ ) as opposed to poly log b-a γ , where log * t is the number of times log has to be applied to t to get it to be ≤ 1). 5 To prove Theorem E.15, we first relate the Wasserstein distance of interest (between the true distribution P and the algorithm's output distribution P DP to a quantity related to an appropriately chosen restriction. Let q α represent the α-quantile of P and qα represent the α-quantile of Pn and qα represent the α-quantiles of P DP . We also note that all these distributions (and others that will come up in the proof) are bounded distributions over the real line and so we can freely apply the triangle inequality for Wasserstein distance, and the cumulative distribution formula for Wasserstein distance (Lemma A.3). The proof of the main theorem will follow from the following lemmas (all proved in Appendix J).Lemma E.16. Let C ′′ > 0 be a sufficiently large constant, and let n > 0 be sufficiently large. With probability at least 1 -O(β) over the randomness in data samples and Algorithm 5,W(P, P DP ) ≤ W(P | q 1 k ,q 1-1 k , P DP | q 1 k ,q 1-1 k ) + C ′′ W(P, P | q 1 k ,q 1-1 k ).Lemma E.17 (Wasserstein in terms of quantiles). For all datasets x (with data entries in [a, b]), with probability at least 1 -β over the randomness of Algorithm 5, we have thatW( Pn | q 1 k ,q 1-1 k , P DP | q 1 k ,q 1-1 k ) ≤ 2 k q 1-1/k -q 1/k ,where Pn is the uniform distribution over x. Now, we argue about the concentration of the Wasserstein distance between restrictions of the empirical distribution and restrictions of the true distribution.Claim E.18. Fix β ∈ (0, 1) and sufficiently large constants c 3 , c 6 . Let n > 0 be sufficiently large such that n > log n/β (as in Theorem E.15). For all k such that 1 k > c 3 log n β n , with probability at least 1 -O(β) over the randomness in the data,W(P | q 1 k ,q 1-1 k , Pn | q 1 k ,q 1-1 k ) ≤ c 6 log n β • E[W(P | q 1 k ,q 1-1 k , Pn | q 1 k ,q 1-1 k )].Now, we give the proof of our main theorem. 5 The theorem would be of the same form as Theorem E.15, except that Algorithm 5 would be (ε, δ)-DP, Theorem E.15. Using Lemma E.16, Claim E.18 and the triangle inequality, we have that with probability at least 1 -O(β) over the randomness of the data and the algorithm,withW(P, P DP ) ≤ W(P | q 1 k ,q 1-1 k , P DP | q 1 k ,q 1-1 k ) + C ′′ W(P, P | q 1 k ,q 1-1 k ) ≤ W( Pn | q 1 k ,q 1-1 k , P DP | q 1 k ,q 1-1 k ) + W( Pn | q 1 k ,q 1-1 k , P | q 1 k ,q 1-1 k ) + C ′′ W(P, P | q 1 k ,q 1-1 k ) ≤ W( Pn | q 1 k ,q 1-1 k , P DP | q 1 k ,q 1-1 k ) + c 6 log n β E W( Pn | q 1 k ,q 1-1 k , P | q 1 k ,q 1-1 k ) + C ′′ W(P, P | q 1 k ,q 1-1 k )Finally, applying Lemma E.17 and taking a union bound over failure probabilities, we get that with probability at least 1 -O(β) over the randomness of the data and the algorithm,W(P, P DP ) ≤ 2 k q 1-1/k -q 1/k + c 6 log n β E W( Pn | q 1 k ,q 1-1 k , P | q 1 k ,q 1-1 k ) + C ′′ W(P, P | q 1 k ,q 1-1 k )as required.Below we describe the experiment referenced in the introduction.The distribution: We have taken a distribution on [0, 999], which is concentrated on two points 430 and 440, with p 430 = 1 3 and p 440 = 2 3 . These algorithms have been run with n = 1600 samples from this distribution.Minimax Optimal Algorithm: The minimax-optimal algorithm here is the algorithm PSMM from [HVZ23] that considers a fixed partitioning of the interval into Ω(m 1 d ) equal intervals and places the empirical mass in each interval on an arbitrary point in each interval. Here we consider this algorithm with ε = ∞, so that no noise is added. We have run it here with K = 40 buckets.Instance-optimal Algorithm: The instance-optimal algorithm finds k quantiles as in Algorithm 5. In this particular implementation, we used the recursive exponential mechanism of [KSS22] , but we expect other quantile algorithms would work similarly. In this particular case, we use k = 10 quantiles with ε = 1.G Proofs for Section D Theorem D.6. For every level ℓ ∈ [D T ], define the neighborhood of P ℓ asN ℓ : ∆([N ℓ ]) → P(∆([N ℓ ])) by N ℓ (P ℓ ) = {Q ℓ | D ∞ (P ℓ , Q ℓ ) ≤ ln 2}. Then, R N ,n,ε (P ) ≥ max ℓ∈[D T ] r ℓ • R N ℓ ,n,ε (P ℓ ),where the error of P is measured in the Wasserstein distance and P ℓ is measured in the TV distance.Proof of Theorem D.6. Given a distribution P , letA * P = arg min A is ε-DP max Q∈N (P ) E D∼P n [W(P, A(D))] so R N ,n,ε (P ) = max Q∈N (P ) E D∼Q n [W(P, A * P (D))]. Let ℓ ∈ [D T ].We want to define an algorithm A * P ℓ on the distributions in N ℓ (P ℓ ) that achieves maximum error rate 1 r ℓ R N ,n,ε (P ). Define a randomised function g P which given a node ν ℓ at level ℓ, g P (ν ℓ ) is sampled from the distribution P restricted to the leaf nodes that are children of ν ℓ . Given a set of nodes at level ℓ, define g P (D) to be the set where g D is applied to each set element individually. Then define A * P ℓ (D) = (A * P (g P (D))) ℓ . Since g P is applied individually to each element in D, A * P ℓ is ε-DP. Given a distribution Q ℓ ∈ N ℓ (P ℓ ), define a distribution Q on the leaves of the tree as follows:Q(ν) = Q ℓ (ν ℓ ) P ℓ (ν ℓ ) * P (ν),where ν ℓ is the parent node of ν at level ℓ. Note Q ∈ N (P ),g P (Q ℓ ) = Q and Q ℓ = Q ℓ . Now, T V (Q ℓ , A * P ℓ (D)) = T V (Q ℓ , (A * P (g P (D)) ℓ ) ≤ 1 r ℓ ℓ ′ ∈[D T ] r ℓ ′ T V (Q ℓ ′ , (A * P (g P (D)) ℓ ′ ) = 1 r ℓ W(Q, A * P (g P (D)))where the first inequality follows by definition of A * P ℓ and the fact Q ℓ = Q ℓ . Since g P (Q ℓ ) = Q, this implies that for all distributions in N ℓ (P ℓ ),E D∼Q ℓ T V (Q ℓ , A * P ℓ (D)) ≤ E D∼Q 1 r ℓ W(Q, A * P (D)) ≤ 1 r ℓ R N ,n,ε (P ),which implies for all levels ℓ, R N ℓ ,n,ε (P ℓ ) ≤ 1 r ℓ R N ,n,ε (P ) and so we are done.Lemma G.1. [A extension of (ε, δ)-DP Assouad's method [ASZ21] ] Let k 0 , k 1 , • • • be a sequence of natural numbers such that s k s < ∞, ε > 0 and δ ∈ [0, 1]. Given a family of distributions P ⊂ ∆(X ) on a space X , a parameter θ : P → M where M is a metric space with metric d, suppose that there exists a set V ⊂ P of distributions indexed by the product of hypercubesEQUATIONFor each coordinate s ∈ N, j ∈ [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) = 2 |E k0 × E k1 × • • • | u∈E k 0 ×E k 1 ×•••:u s j =+1 p u , p -(s,j) = 2 |E k0 × E k1 × • • • | u∈E k 0 ×E k 1 ×•••:u s j =-1 p u ,and let ϕ s,j : X n → {-1, +1} be a binary classifier. Thenmin A is (ε,δ)-DP max p∈V R A,n (p) ≥ 1 2 s τ s ks j=1 min ϕs,j is (ε,δ)-DP ( Pr X∼p n +(s,j) (ϕ s,j (X) ̸ = 1)+ Pr X∼p n -(s,j) (ϕ s,j (X) ̸ = -1)),where the min on the LHS is over all (ε, δ)-DP mechanisms, and on the right hand side is over all (ε, δ)-DP binary classifiers. Moreover, if for all s ∈ N, j ∈ [k s ], there exists a coupling (X, Y ) between p n +(s,j) andp n -(s,j) with E[d Ham (X, Y )] ≤ D s , then min A is (ε,δ)-DP max p∈V R A,n (p) ≥ s k s τ s 2 (0.9e -10εDs -10D s δ)Proof of Lemma D.8. We will follow the proof of Theorem 3 in [ASZ21] . Given an estimator A, define a classifier A * by projecting on the product of hypercubes soA * (X) = arg min u∈(E k 0 ×E k 1 ×••• ) d(A(X), θ(p u )).By the triangle inequality and the definition of A * , for any p ∈ V,d(θ(p A * (X) ), θ(p)) ≤ d(A(X), θ(p A * (X) )) + d(A(X), θ(p)) ≤ 2d(A(X), θ(p)).Therefore, we can restrict to a lower bound on the performance of DP classifiers:EQUATIONin which case P ℓ (nu) ≤ 4 3P ℓ (ν) ln(n/β) n , as required. Similarly, if P ℓ (ν) > 1/2 then truncation will only occur if 1 -P ℓ (ν) ≤ 4 3(1-P ℓ (ν)) ln(n/β) n , as required.We will now show that these conditions hold simultaneously with probability at least 1 -3β for all the nodes at level ℓ. If P ℓ (ν) ≤ 1 en then using the multiplicative form of Chernoff bound,Pr ( Pℓ ) n (ν) ≥ 3 ln(n/β) n = Pr ( Pℓ ) n (ν) ≥ 1 + 3 ln(n/β) P ℓ (ν)n -1 P ℓ (ν) ≤   e 3 ln(n/β) nP ℓ (ν) -1 ( 3 ln(n/β) nP ℓ (ν) ) 3 ln(n/β) nP ℓ (ν)   P ℓ (ν)n ≤ enP ℓ (ν) 3 ln(n/β) 3 ln(n/β) ≤ P ℓ (ν)n e 3 ln(n/β) 3 ln(n/β) (nP ℓ (ν)) 3 ln(n/β)-1Firstly, since ln(n/β) ≥ 1, e 3 ln(n/β) 3 ln(n/β) ≤ 1. Further, nP ℓ (ν) ≤ 1/e and 3 ln(n/β) -1 ≥ ln(n/β) so (nP ℓ (ν)) 3 ln(n/β)-1 ≤ (1/e) ln(n/β) = β/n. Therefore,EQUATIONLet S = {x ∈ [N ℓ ] | P ℓ (x) < 1/(en)} then using a union bound and Eqn (12) we havePr ∃x ∈ S s.t. (P ℓ ) n (x) ≥ 2 √ 2 log(n/β) n ≤ x∈S P ℓ (ν)β ≤ βThere exist at most n elements in [N ℓ ] that do not belong in S. We will prove that, independently, each of these elements satisfy the required condition with probability ≤ 2β/n then a union bound proves the final result. If P ℓ (ν) ∈ [ 3 ln(n/β) n , 1 -3 ln(n/β) n ] then using the multiplicative form of Chernoff bound (If X i are all i.i.d. and 0 < δ < 1, then Pr(|n i=1 X i -nE[X 1 ]| ≥ δnE[X 1 ]) ≤ 2e -δ 2 nE[X1]/3 ), Pr | (P ℓ ) n (x) -P ℓ (x)| ≥ 3P ℓ (x) log(n/β) n = Pr | (P ℓ ) n (x) -P ℓ (x)| ≥ 3 log(n/β) P ℓ (x)n P ℓ (x) ≤ 2e - 3 log(n/β) P ℓ (x)n P ℓ (x)n 3 = 2β/n.Next, if P ℓ (ν) ≤ 3 ln(n/β) n then using the additive form of Chernoff bound (If X i are all i.i.d. and ε ≥ 0, then Pr EQUATIONG P (ν) + Lap( 1 εn ) ≥ G P (ν) -4 3 G P (ν) log(n/β) n - ln(2/β) εn ≥ 1 2 G P (ν) - log(2/β) εn ≥ 1 εn + log(2/β) nEach level has at most 2εn in γ P 1 2εn so we query at most 4εn nodes in the tree when running LocateActiveNodes since each node has at most 2 children. Therefore, we can set X = 4εnD )as required.Lemma G.7. For any real-valued function G on the nodes of the HST such that G(ν 0 ) = 1 where ν 0 is the root node and given any distribution P ,W(P, Projection(G)) ≤ 4W(G P , G).Proof of Lemma D.17. We first note that for any pair of sequences of real valuesa 1 , • • • , a k and b 1 , • • • , b k , and constant A such that i a i ̸ = 0, | A a i a i -b i | ≤ | A a i a i -a i |+|a i -b i | = |A- a i |+ |a i -b i | ≤ |A- b i |+2 |a i -b i |. Also if a i = 0 then | A k -b i | ≤ | A k -i b i k |+| i b i k -b i | = |A- b i |+2 |b i | = |A- b i |+2 |a i -b i |Let Ḡℓ be the function Ḡ after only levels 0, • • • , ℓ have been updated. So Ḡℓ matches Ḡℓ-1 on all levels except ℓ. Let ν be a node in the ℓth level of the HST. If we suppose the sum is over the normalised children of a node ν, A = Ḡℓ-1 (ν), and for all the children ν ′ of ν, a i = G(ν ′ ) and b i = G P (ν ′ ), we can see that the contribution to the Wasserstein distance by the children increases by an additive factor of |G ℓ-1 (ν) -G P (ν)|. Iterating, we can see thatW(P, Projection(G)) ≤ 2 D T ℓ=0 ν at level ℓ (r ℓ +r ℓ+1 • • • r D T )|G(ν)-G P (ν)| ≤ 4 D T ℓ=0 ν at level ℓ r ℓ |G(ν)-G P (ν)|,which is 4 times the wasserstein distance.Theorem H.1. Given any ε > 0, and a distribution P , and let n ′ = 5 4 min{W ( 0.45ε δ ),0.6} n, then for all (ε, δ)-DP algorithms A ′ , there exists a distribution Q ∈ N (P ) such that with probability 1 -(D T log n + 4D T εn)β,W(Q, Q ε,n ′ ) ≤ Õ(E X∼Q n ,A ′ (W(A ′ (X), Q))),where Q ε,n ′ is the output of PrivDensityEstTree(Q) with n ′ samples.Proof. First, let us obtain a slightly simpler upper bound on W(P, Pε ). From eqn (9) in the proof of Theorem D.13 we have that for each level ℓ, W(P ℓ , ( Pε ) ℓ ) ≤ 2 W((G P ) ℓ , ( G P ) ℓ ) + W(( G P ) ℓ , ( G P |γ ε ) ℓ ) + W(( G P |γ ε ) ℓ , ( G Pn,γε ) ℓ ) ,≤ 2C ln(n/β)    ν min G P (ν)(1 -G P (ν)), G P (ν)(1 -G P (ν)) n + ν / ∈γ P ( 1 2εn ) G P (ν) + (|γ P 1 2εn | -1) 1 εn    ,Further, by Theorem D.7 and Theorem D.6, given ε > 0 and δ ∈ [0, 1], let κ = 1 10εn min{W 0.45ε δ , 0.6} where W (x) is the Lambert W function so W (x)e W (x) = x. Given a distribution P , there exists a constant C ′ such thatR N ,n,ε (P ) ≥ C ′ D T   ν min G P (ν)(1 -G P (ν)), G P (ν)(1 -G P (ν)) n + ν / ∈γ P (2κ) G P (ν) + (|γ P (2κ) | -1)κ   Let Q ∈ N (P ), then γ P 1 εn ⊂ γ Q 1 2εn ⊂ γ P 1 4εn so ν / ∈γ Q ( 1 2εn ) G Q (ν) + (|γ Q 1 2εn | -1) 1 εn = ν / ∈γ P ( 1 4εn ) G Q (ν) + ν∈γ P ( 1 4εn )\γQ( 1 2εn ) G Q (ν) + (|γ Q 1 2εn | -1) 1 εn ≤ ν / ∈γ P ( 1 4εn ) 2G P (ν) + ν∈γ P ( 1 4εn )\γQ( 1 2εn ) 1 εn + (|γ Q 1 2εn | -1) 1 εn ≤ ν / ∈γ P ( 1 4εn ) 2G P (ν) + (|γ P 1 4εn | -1) 1 εn . Now, let n ′ = 5 4 min{W ( 0.45ε δ ),0.6} n ≥ n so for all Q ∈ N (P ), W(Q, Q ε,n ′ ) ≤ Õ    ν min G Q (ν)(1 -G Q (ν)), G Q (ν)(1 -G Q (ν)) n ′ + ν / ∈γ Q ( 1 2εn ′ ) G Q (ν) + (|γ Q 1 2εn ′ | -1) 1 εn ′    ≤ Õ    ν 2 min G P (ν)(1 -G P (ν)), G P (ν)(1 -G P (ν)) n ′ + ν / ∈γ P ( 1 4εn ′ ) 2G P (ν) + (|γ P 1 4εn ′ | -1) 1 εn ′    = Õ   ν 2 min G P (ν)(1 -G P (ν)), G P (ν)(1 -G P (ν)) n + ν / ∈γ P (2κ) 2G P (ν) + (|γ P (2κ) | -1) 1 εn   ≤ Õ min A ′ max Q ′ ∈N (P ) E X∼Q ′n ,A ′ (W(A ′ (X), Q ′ )) .As in Proposition B.5, since N (P ) is compact, for all A ′ , there exists a specificQ * ∈ N (P ) such that W(Q * , Q * ε,n ′ ) ≤ Õ(E X∼(Q * ) n ,A ′ (W(A ′ (X), Q * )))Estimating appropriately chosen quantiles is the main part of our algorithm for approximating the distribution over R in Wasserstein distance, and so in this section, we describe some known differentially private algorithms for this task and derive some corollaries that we use extensively in our application. We will use F to represent CDF functions, with F P representing the CDF of distribution P . We start by stating an important theorem on private CDF estimation. This follows from a use of the binary tree mechanism [CSS11, DNPR10] . A version of this theorem for approximate differential privacy is described in a survey by Kamath and Ullman [KU20, Theorem 4.1]. The version presented here for pure differential privacy follows from a very similar argument, except using Laplace Noise instead of Gaussian noise (and basic composition instead of advanced composition to analyze privacy). Their accuracy was also in expectation, but a similar analysis yields a high probability bound, as in the theorem below.Theorem I.1. [KU20, Theorem 4.1] Let ε, β ∈ (0, 1], let D be an ordered, finite domain, and let x ∈ D n be a dataset. Let Pn be the uniform distribution on x. Then, there exists an ε-DP algorithm A CDF that on input x and the domain D outputs a vector G over D such that with probability at least 1 -β over the randomness of A CDF :∥G -F Pn ∥ ∞ = O log 3 |D| β εn .CDF estimation is intimately related to quantile estimation, and we use the following quantitative statement that will follow from a simple application of Theorem I.1.Theorem I.2. Fix any n > 0, ε, β ∈ (0, 1], a, b ∈ R, and γ < b -a ∈ R such that b-a γ is an integer. Let C be a sufficiently large constant. Then, there exists an ε-DP algorithm A quant , that on input interval end points a, b, granularity γ, x = (x 1 , . . . , x n ) ∈ {a, a + γ, . . . , b -γ, b} n , and desired quantile values α ∈ (0, 1) k , outputs quantiles q ∈ {a, a + γ, . . . , b -γ, b} k such that with probability at least 1 -β over the randomness of A quant , for all r ∈EQUATIONPr y∼ Pn (y < qr ) < α r + C log 3 b-a βγ εn ,where Pn is the uniform distribution on the entries of x.Proof. Algorithm A quant operates by running the algorithm A CDF referenced in Theorem I.1 on x and domain {a, a + γ, . . . , b -γ, b}, and postprocessing its outputs to get quantile estimates as follows. For every quantile α r that we are asked to estimate, A quant simply scans the vector G output by algorithm A CDF in order, and outputs the first domain element whose CDF estimate in G crosses α r . Conditioned on the accuracy of the CDF estimation algorithm G, we have that this output qr satisfiesα r -F Pn (q r ) ≤ C log 3 b-a βγ εn .Additionally, since qr is the first domain element whose estimate in G crosses α r , we also have that P r y∼ Pn (y < qr ) < α r + C log 3 b-a βγ εn .Hence, with probability at least 1 -β, we have this property for all r ∈ [k].We now state a corollary of this theorem that we will use extensively in our presentation.Corollary I.3. Fix any ε, β ∈ (0, 1], a, b ∈ R, and γ < b -a ∈ R such that b-a γ is an integer. Let n ∈ N > 4c2 log 4 ( b-a βγε ) ε , such that k set to ⌈ εn 4c3log 3 b-a βγ log n β⌉ is an integer greater than or equal to 1, where c 2 and c 3 are sufficiently large constants.6 Then, there exists an ε-DP algorithm A quant (the same one referenced in Theorem I.2), that on input interval end points a, b, granularity γ, x = (x 1 , . . . , x n ) ∈ {a, a + γ, . . . , b -γ, b} n , and desired quantile values α = {1/2k, 3/2k, 5/2k, . . . , (2k -1)/2k}, outputs quantiles q ∈ {a, a + γ, . . . , bγ, b} k such that with probability at least 1 -β, for all r ∈ [k],q 2r-1 2k -1 4k ≤ qr ≤ q 2r-1 2k + 1 4k ,where Pn is the uniform distribution on the entries of x and for all p ∈ (0, 1), qp is the p-quantile of Pn .Proof. First, note that k is set such that 1 4k ≥ C This implies that for all r ∈ [k], qr ≤ q 2r-1 2k + 1 4k .J Proofs in Section E J.1 Omitted Proofs in Section E.1.1Proof of Lemma E.6. We evaluate the various terms in Theorem E.5.We start by evaluating τ (P,Q) = max{ R max(f P (t) -e ε f Q (t), 0)dt, R max(f Q (t) - e ε f P (t), 0)dt }.Consider the first term in the outer maximum. For all t ∈ [L(P ), q 1/k ), we have that f Q (t) = 1 2 f P (t). For all other t, one can see that the value of the integrand is 0. Hence, the value of the first term isq 1/k L(P ) max(f P (t) -e ε 2 f P (t), 0)dt = max{ 1 -e ε 2 1 k , 0} ≤ 1 2k. Now, consider the second term in the outer maximum. For all t < q 1-1 k , the value of the integrand is 0. For all q 1-1 k ≤ t ≤ q 1 , the value of the integrand is max{ 3 2 -e ε f P (t), 0}. Hence, the second term is max{ 3 2 -e ε 1 k , 0} ≤ 1 2k . Put together, we get that τ (P, Q) ≤ 1 2k . When ε ≥ ln 2, we have that 1 -e ε 2 ≤ 0, and so we have that the largest value of ε′ ∈ [0, ε] that makes R max(f Q (t) -e ε ′ f P (t), 0)dt = τ (P, Q) = 0, is ε ′ = ε. When ε < ln 2, we have that the value of ε ′ that makes R max(f Q (t) -e ε ′ f P (t), 0)dt = max{ 3 2 -e ε 1 k , 0} = 1 -e ε 2 1 k , is ε ′ = ln 1+e ε 2 .Finally, we describe the distributions P ′ and Q ′ and compute the squared Hellinger distance between them. There are two cases, based on the range of ε. First, consider ε ≥ ln 2. First, we calculate P ≡ min{e ε Q, P }. This value is equal to min{e ε /2, 1}f P (t) = f P (t) for t < q 1 k (P ), and is also equal to f P (t) for q 1 k ≤ t ≤ q 1 . Similarly, consider Q ≡ min{e ε ′ P, Q} = min{e ε P, Q}; it is equal to f P (t) 2 for t < q 1 k (P ), and is equal to f P (t) for q 1 k ≤ t ≤ q 1-1 k . It is also equal to min(e ε , 32 )f P (t) = 3 2 f P (t) for q 1 k ≤ t ≤ q 1 . Since τ (P, Q) = 0, and by the above calculations, we have that P ′ = P , and Q ′ = Q. Upper bounding the squared Hellinger distance between P ′ and Q ′ by the TV distance (See Lemma A.7), we get that H2 (P ′ , Q ′ ) = H 2 (P, Q) ≤ T V (P, Q) = 1 2k ≤ ε 2(ln 2)k (where we have used that ε ≥ ln 2).Next, consider ε < ln 2. First, consider P ≡ min{e ε Q, P }. This value is equal to min{e ε /2, 1}f P (t) = e ε 2 f P (t) for t < q 1 k (P ), and is also equal to f P (t) for q 1 k ≤ t ≤ q 1 . Similarly, consider Q ≡ min{e ε ′ P, Q} = min{ 1+e ε 2 P, Q}; it is equal to 1 2 f P (t) for t < q 1 k (P ), and is equal to fP (t) at q 1 k ≤ t ≤ q 1-1 k . It is also equal to min{ 1+e ε 2 , 3 2 }f P (t) = 1+e ε 2 f P (t) at q 1-1 k ≤ t ≤ q 1 . Note that τ (P, Q) = 1 -e ε 2 1k . P ′ and Q ′ are the distributions created by normalizing P and Q by dividing by a factor of 1 -τ (P, Q). Now, we upper bound the squared Hellinger distance between P ′ and Q ′ by the TV distance (See Lemma A.7), to get thatH 2 (P ′ , Q ′ ) ≤ T V (P ′ , Q ′ ) = O( ε k). Substituting into the lower bound for sample complexity of distinguishing P and Q, this tells us that for all ε ∈ (0, 1], SC ε (P, Q) = Ω 1 ε• 1 k = Ω(k/ε).Proof of Lemma E.7. Note that P has bounded expectation (and hence, so does Q). Hence, we can use the following form of the Wasserstein distance:W(P, Q) = R |F P (t) -F Q (t)|dt.where the last inequality is by pulling the first two terms from the summation in second term to the summation in the first term, and using the fact that for j = log n -2, j = log n -1, we have that1 4•2 j ≥ 1 2 √ 2 1 √ 2 j nProof of Lemma E.12. We first state a theorem of Bobkov and Ledoux [BL19] .Theorem J.1 (Theorem 3.5, [BL19] ). There is an absolute constant c > 0, such that for all distributions P over R, for every n ≥ 1,c(A n + B n ) ≤ E[W(P, Pn ] ≤ A n + B n .whereEQUATIONB n = 1 √ n F (t)[1-F (t)]≥ 1 4n F (t)[1 -F (t)]dt.Now, we are ready to prove the main theorem. Fix natural number i ≥ 2. Restricted to t ≤ q 1/2 , F P (t)(1 -F P (t)) is an increasing function, and hence for t ∈ [q 1/2 i , q 1/2 i-1 ], we have thatF P (t)(1 -F P (t)) ≤ 1 2 i-1 [1 -1 2 i-1 ].Similarly, restricted to t > q 1/2 , F P (t)(1 -F P (t)) is a decreasing function, and hence for t ∈[1 -q 1/2 i-1 , q 1-1/2 i ], we have that F P (t)(1 -F P (t)) ≤ 1 2 i-1 [1 -1 2 i-1 ].Using this, we can now upper bound the expected Wasserstein distance between P and its empirical distribution using Theorem J.1. Hence, we upper bound the terms B n and A n . We start by upper bounding B n . Note that for all t ̸ ∈ [q 1 4n , q 1-1 4n ], we have that F P (t)(1 -F P (t)) ≤ 1 4n . Hence,B n = 1 √ n F P (t)[1-F P (t)]≥ 1 4n F P (t)[1 -F P (t)]dt ≤ 1 √ n q 1-1 4n q 1 4n F P (t)[1 -F P (t)]dt ≤ log 4n i=2 1 √ n q 1/2 i-1 q 1/2 i F P (t)[1 -F P (t)]dt + q 1-1/2 i q 1-1/2 i-1 F P (t)[1 -F P (t)]dt ≤ log 4n i=2 1 √ n q 1/2 i-1 q 1/2 i 1 2 i-1 1 - 1 2 i-1 dt + q 1-1/2 i q 1-1/2 i-1 1 2 i-1 1 - 1 2 i-1 dt = log 4n i=2 1 √ n 1 2 i-1 1 - 1 2 i-1 q 1/2 i-1 -q 1/2 i + q 1-1/2 i -q 1-1/2 i-1 ≤ log 4n i=2 2 √ 2 i n q 1-1/2 i -q 1/2 i = log n-1 i=2 2 √ 2 i n q 1-1/2 i -q 1/2 i + log 4n i=log n 2 √ 2 i n q 1-1/2 i -q 1/2 i ≤ log n-1 i=2 2 √ 2 i n q 1-1/2 i -q 1/2 i + log 4n i=log n 4 2 i q 1-1/2 i -q 1/2 i ,where the last inequality is because for i ≤ log(4n), we have that 1 n ≤ 4 2 i . Next, we bound A n . Note that for all t ≥ q 1/2n and for all t ≤ q 1-1 2n , we have that F P (t)(1 -F P (t)) ̸ ≤ 1 4n . Hence,A n = 2 F P (t)[1-F P (t)]≤ 1 4n F P (t)[1 -F P (t)]dt ≤ 2   q 1 2n -∞ F P (t)[1 -F P (t)]dt + ∞ q 1-1 2n F P (t)[1 -F P (t)]dt   = ∞ i=1+log 2n 2 q 1/2 i-1 q 1/2 i F P (t)[1 -F P (t)]dt + q 1-1/2 i q 1-1/2 i-1 F P (t)[1 -F P (t)]dt ≤ ∞ i=1+log 2n 2 q 1/2 i-1 q 1/2 i 1 2 i-1 1 - 1 2 i-1 dt + q 1-1/2 i q 1-1/2 i-1 1 2 i-1 1 - 1 2 i-1 dt = ∞ i=1+log 2n 2 2 i-1 1 - 1 2 i-1 q 1/2 i-1 -q 1/2 i + q 1-1/2 i -q 1-1/2 i-1 ≤ ∞ i=1+log 2n 4 2 i q 1-1/2 i -q 1/2 iThen, using the upper bound in Theorem J.1, substituting in the bounds for A n and B n , and simplifying, we get the claim.Proof of Claim E.13. By the definition of Wasserstein distance and restrictions of distributions, we have that W( Pn | q 1 k ,q 1-1 k , P | q 1 k ,q 1-1 k ) = b a F Pn|q 1 k ,q 1-1 k (t) -F P |q 1 k ,q 1-1 k (t) dt = q 1-1 k q 1 k F Pn|q 1 k ,q 1-1 k (t) -F P |q 1 k ,q 1-1 k (t) dt = q 1-1 k q 1 k F Pn (tNote that for all t ∈ [q 1/k , q 1-1/k ], we have that the cumulative distribution functions of P and its restricted version are identical and likewise for P DP . Additionally, the cumulative density functions for the restricted versions of the two distributions are identical to each other outside of this interval. Hence, we can simplify the middle term in the RHS of the inequality above as follows:q 1-1/k (P )t=q 1/k |F P (t) -F P DP (t)|dt = W(P | q 1 k ,q 1-1 k , P DP | q 1 k ,q 1-1 k )Next, we reason about the remaining terms.Consider the term q 1/k t=a |F P (t) -F P DP (t)|dt. First, condition on the event in Theorem I.3 (on the accuracy of the private quantiles for the empirical distribution), which tells us that with probability at least 1 -β, we have for all r ∈ [k], thatEQUATIONwhich implies in particular that q1/4k ≤ q1/2k ≤ q3/4k .Next, we argue that q1/4k ≥ q 1/8k with high probability. By the definition of quantiles, we have that P r y∼P (y < q 1/8k ) < 1 8k . The number of entries in the dataset x less than q 1/8k is hence a Binomial with mean less than n 8k , and hence, we have by Theorem J.2 (with δ set to 0.9) that with probability at least 1 -β, the number of entries in the dataset less than q 1/8k is at most 1.9 n 8k < n 4k , which means the total mass less than q 1 8k in the empirical distribution is less than 1 4k . This implies that q1/4k ≥ q 1/8k by the definition of quantiles.Additionally, note that for all t < q 1/k , F P (t) < 1 k . The number of entries in the dataset x that are less than q 1/k is hence a Binomial with success probability less than 1 k . By Theorem J.2, we can again argue that with probability at least 1 -β, there is a constant c ′ such that the total mass of the empirical distribution on values less than q 1/k is less than c ′ k . Hence, q 1/k ≤ qc ′ /k . This implies by Equation 15, that q 1/k ≤ qc/k for some constant c. Hence, for all t < q 1/k , we have that F P DP (t) ≤ c k . Hence, taking a union bound, with probability at least 1 -O(β), ).q 1/Taking a union bound to ensure that all terms in Equation 14 are bounded as required, the proof is complete.Proof of Lemma E.17. First, we condition on the event in Corollary I.3 (on the accuracy of differentially private quantile estimates) that for all r ∈ [k], q 2r-1 2k -1 4k ≤ qr ≤ q 2r-1 2k + 1 4k , note that this event happens with probability at least 1 -β over the randomness of the algorithm.Observe that this implies that F DP increases by 1 k somewhere in the range [q 2r-1 2k -1 4k , q 2r-1 2k + 1 4k ] (for all r ∈ [k]) and remains constant outside these intervals. Now, we show that for all t ∈ [a, b], we have that |F P DP (t) -F Pn (t)| ≤ 2 k . If there exists t ∈ [a, q 1 4k ), we have that F P DP (t) = 0, and F Pn (t) ≤ 1 4k , which implies that |F DP (t)-F Pn (t)| ≤ 1 4k . If there exists no such t, then we have that a = q 1 4k , and the corresponding interval collapses to a single point (which will fall in another interval considered below).Next, fix any r ∈ [k -1]. Note that if there exists t ∈ [q 2r-1 2k -1 4k , q 2r+1 2k -1 4k ), we have for all such t that r-1 k ≤ F DP (t) < r k , and 2r-1 2k -1 4k ≤ F Pn (t) ≤ 2r+1 2k + 1 4k . This implies that for all such t, |F DP (t) -F Pn (t)| ≤ 2 k . If there exists no such t, then we have that q 2r-1 2k -1 4k = q 2r+1 2k -1 4k, and this r is not relevant since the corresponding interval collapses to a single point (that is considered in another interval).Finally, for t ∈ [q 2k-1 2k, b], we have that F P DP (t) ≥ 1 -1 k , and F Pn (t) ≥ 1 -1 2k , so we have that |F DP (t) -F Pn (t)| ≤ 1 k . Note that every t ∈ [a, b] is considered in some interval above and hence we have shown that for all t ∈ [a, b], we have that |F P DP (t) -F Pn (t)| ≤ 2 k . Finally, using the formula for Wasserstein distance (and the definition of a restriction), we have thatEQUATION≤ q 1-1 k q 1 k 2 k dt (18) ≤ 2 k q 1-1/k -q 1/k (19)Before the proof of Claim E.18, we state the following variance-dependent version of the DKW inequality that uniformly bounds the absolute difference in CDFs between the true and empirical distribution. Theorem J.3 (See for example Theorem 1.2 in [BM23] ). Fix n > 0. There are absolute constants c 0 , c 1 such that for all ∆ ≥ c0 log log n n , Pr sup t:F P (t)(1-F P (t))≥∆F P (t) -F Pn (t) ≥ ∆ • F (t)(1 -F (t) ≤ 2e -c1∆nWe also state the following lemma on Binomial random variables, which is a simple consequence of a Lemma by Bobkov and Ledoux [BL19] . Lemma J.4 (Lemma 3.8 in [BL19] ). Let S n = n i=1 η i be the sum of n independent Bernoulli random variables with Pr[η i = 1] = p and Pr[η i = 0] = q = 1 -p (for all i). Also assume p ∈ [ 1 n , 1 -1 n ]. Then, for some sufficiently small constant c, c√ npq ≤ E[|S n -np|] ≤ √ npqProof of Claim E.18. Now, by the formula for Wasserstein distance, the definition of restriction, and Fubini's theorem, we have thatE[W(P | q 1 k ,q 1-1 k , Pn | q 1 k ,q 1-1 k )] = E q 1-1 k q 1 k F P (t) -F Pn (t) dt = q 1-1 k q 1 k E F P (t) -F Pn (t) dtBy Lemma J.4, using the fact that F Pn (t) = n i=1 1[x i ≤ t], where each term in the sum is an independent Bernoulli random variable with expectation F P (t), with q 1 k ≤ t < q 1-1 k (ensuring that the conditions of the lemma are met), we get that E F P (t) -F Pn (t) ≥ c F P (t)[1-F P (t)] n , which givesE[W(P | q 1 k ,q 1-1 k , Pn | q 1 k ,q 1-1 k )] ≥ c q 1-1 k q 1 k F P (t)[1 -F P (t)] n dtNow, consider the random variable W(P | q 1 k ,q 1-1k , Pn | q 1 k ,q 1-1 k ). Note that 1 k ≥ c3 log n β n(for an appropriately chosen c 3 ), and so we are in the regime where we can apply Theorem J.3 for an appropriately chosen ∆.In particular, we have that for t ∈ [q 1 k , q 1-1 k ),F P (t) ∈ [ 1 k , 1 -1 k ).Setting ∆ = , and ∆ ≤ 1 2k (the second inequality for sufficiently large c 3 ). In particular, this implies for t ∈ [q 1 k , q 1-1 k ), F P (t) ∈ [2∆, 1 -2∆), which implies that F P (t)(1 -F P (t)) ≥ ∆, as long as n > c 4 log n β for some sufficiently large constant c 4 . Now, using Theorem J.3, we have that with probability at least 1 -2e -c1log n β c 1 n n ≥ 1 -O(β), sup t∈[q 1 k ,q 1-1 k ) F P (t) -F Pn (t) ≤ log n β c 1 n F P (t)(1 -F P (t))Condition on this for the rest of the proof. Then, we can write the following set of equations.W(P | q 1 k ,q 1-1 k , Pn | q 1 k ,q 1-1 k ) = q 1-1 k q 1 k |F P (t) -F Pn (t)|dt ≤ q 1-1 k q 1 k log n β c 1 n F P (t)(1 -F P (t))dt ≤ c 5 log n β q 1-1 k q 1 k F P (t)(1 -F P (t)) n dt ≤ c 6 log n β E[W(P | q 1 k ,q 1-1 k , Pn | q 1 k ,q 1-1 k )]as required.In this subsection, we argue that the instance-optimal algorithm discussed in Section E.2 is also locally-minimal (See Appendix B.2 for a discussion of local minimality).First, we state a corollary of our upper bound for continuous distributions, Theorem E.14. This corollary follows by discretizing the distribution and applying the previous upper bound to the discretized distribution. The parameters of the discretized distribution are related to that of the original distribution via simple coupling arguments. for some sufficiently large constant c 2 . Then, there exists an algorithm, that when given inputs x ∼ P n , 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments? Answer: [No] Justification: [TODO] The experiment is not the focus of our paper (rather, we are focused on theoretical analysis) and is just used to demonstrate one drawback of worst case bounds (which we also discuss theoretically) -hence, we don't comment on statistical significance in detail. Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments? Answer: [Yes] Justification: [TODO]We do not explicitly address this since the experiment can be performed on essentially any laptop since it is simple (both time and memory efficient) and does not require any significant compute. Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 9. Code Of Ethics Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines? Answer: [Yes] Justification: [TODO] We have reviewed the code of ethics in detail-both societal impact guidelines and impact-mitigation measures and are confident that our paper conforms to the code of ethics in every respect.• Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.c.f. related work section for discussion of other beyond worst case analysis approaches for this question and Appendix B for a more in-depth discussion of our approach.√t∈R:f P (t)>0 ( fP (t) -f Q (t)) 2 dt.Note that we use H 2 (P, Q) to represent the squared Hellinger distance. Next, we define total variation distance, which will come up in our high-dimensional results.We 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.The same bounds (and hence all our results in this subsection) can be extended to (ε, δ)-DP (with δ ≤ ε) by using an equivalence of pure and approximate DP for identity and closeness testing [ASZ17, Lemma5].k is set to be sufficiently small in order to relate the accuracy of the quantiles algorithm to a parameter depending on k, and n is set sufficiently large that k is not less than 1. The dependence on β comes up in the proof of Claim E.18. 
 
 Appendix:  


Decision: This paper considers the problem of learning potentially arbitrary distributions with respect to the Wasserstein distance under differential privacy constraints. The authors define a natural notion of instance-optimality in this setting and develop a nearly-instance-optimal learning algorithm (within poly-logarithmic factors). While there was a divergence of scores, after reading the reviews and subsequent discussion, my understanding is that this work is technically and conceptually novel, and clearly meets the acceptance bar at NeurIPS.

Reviews: review: paper_summary: paper_summary: The paper studies distribution estimation under Wasserstein distance while requiring the estimator to be differentially private and "instance-optimal," i.e., being competitive against the algorithm optimal on a distribution neighborhood. The derived class of estimators can achieve comparable performance up to polylogarithmic factors and is extendable to other HST (hierarchically separated tree) metric spaces. Besides the theoretical advancement, the primary conceptual contribution is the new formulation of instance optimality in distribution estimation.
strength_weakness: strength_weakness: - It seems that the reference (optimal) estimator evaluates its performance on both the actual underlying distribution and the worst possible distribution, which differs from the actual by a constant factor of 2 in the distribution density ratio. This formulation significantly weakens the power of the reference distribution since, intuitively, it has to produce something in between that could largely deviate from the actual underlying distribution.
- The poly-log multiplicative factors in the upper bounds do not seem optimal. It would be better if the authors could provide some lower bounds to justify this dependence.
- The formulation is much weaker than some prior works in the discrete distribution setting. For example, OS15 shows that their estimator is comparable to all "natural estimators" on each distribution instance, where the best natural estimator essentially knows the actual distribution (without permutation) and only has to output the same probability estimator for symbols that appear the same number of times (a natural requirement).
questions: questions: Please see the comments in the "weaknesses" section. I appreciate the authors' time and efforts.
limitations: limitations: The authors discussed and compared their results, which seems sufficient given the page limit. Negative societal impact is not applicable here.
paper_summary: paper_summary: The paper focuses on estimating densities while preserving privacy (differential privacy). Preserving privacy comes with trade-offs in terms of accuracy of their estimates, which is measured using the Wasserstein distance in this paper. One of the main conceptual contribution is introducing a notion of instance optimality for this problem. The authors motivate the choice by comparing the notion against some alternatives from related literature. 

Their new definition of instance optimality is similar to the hypothesis testing problem where one has to look at a sequence of samples and decide whether they come from distribution X or Y. Here, instance optimality is defined as the performance of an algorithm that a priori knows that the samples are either from X or Y, and once it successfully determines which, it can return the density. For some real algorithm to be instance optimal, it has to match the performance of this algorithm without the prior knowledge, unto a constant slowdown.  There are some further nuances to the definition that make the definition robust to trivial edge cases.

The paper then considers (DP) estimation with respect to Wasserstein error, which is a well motivated problem in ML. The paper shows improved algorithms for private estimation, that beat the previous baselines wrt to the special case of Discrete distributions and TV metric. 

One other technical insight is in a tighter analysis of existing HST algorithms which are a surprisingly good fit for the problem when considered over R^2.
strength_weakness: strength_weakness: Minor Typos — 
(250) follow

Perhaps not the main focus of the paper, but it would be nice to have all the constants. It encourages implementation, and also constants seem to make or break practical DP algorithms. And further since algorithms are inspired from practical HST algorithms, aren't a some of the source constants already available?
questions: questions: --
limitations: limitations: --
paper_summary: paper_summary: This paper looks at the problem of private instance optimal density estimation in wasserstein distance. The starting point of their research is the recent result that shows that the minimax rate of  eps differentially private density estimation in wasserstein distance scales as (eps*n)^{-1/d}. Therefore, the authors take an instance optimal approach, where given samples from an unknown distribution P the algorithm wishes to achieve a cost of at most alpha times the maximum cost of any other algorithm’s cost in a neighborhood around P, that knows neighborhood. For defining the neighborhood, they consider the class of distributions that are close to P in D infinity distance. Coming up with this notion of instance optimality is claimed to be the main contribution of the paper. They authors provide tight upper and lower bounds for such an instance optimal private density estimation guarantee. The upper bounds mostly use the existing techniques from prior work.

The first result shows that when P is one dimensional and its sample space is [0,1], they can achieve an polygon-factor approximate cost as compared to an algorithm that knows that the unknown distribution is either P or Q, for some distribution Q that exists. Question to the authors, do they mean for every Q or there exists such a Q? If the later is the case, the theorem sounds weak to me.

Similarly, for two dimensional distributions over [0,1]^2, they could get a (log n)^{O(1)} approximation (how is it different from polylog n?) as compared to an algorithm that knows the distribution is ln 2 close to P in D infinity distance.
strength_weakness: strength_weakness: I did not find any theorems that justifies that their polylog n factors are tight or not. I would like to believe the approximation factors are quite pessimistic. Moreover, the results do not apply beyond two dimensions. The choice n’=n/polylog(n) also seems ad-hoc.

I found the results derived in this paper are quite weak for a venue like NeurIPS.
questions: questions: None.
limitations: limitations: None.
paper_summary: paper_summary: The paper addresses the problem of differentially private density estimation using the Wasserstein distance. It approaches this problem in a non-parametric manner, meaning it does not assume any specific distribution. Instead of focusing on worst-case error guarantees, the authors aim to design algorithms that achieve instance-optimal error. In other words, the algorithms provide the best possible error within a small "neighborhood" around any given input.

The authors explore various definitions of "neighborhood" in the context of instance optimality for distribution estimation in both $\mathbb{R}$ and $\mathbb{R}^2$. They present algorithms that achieve performance close to the instance-optimal ones, up to multiplicative polylogarithmic factors. Furthermore, the algorithm for $\mathbb{R}^2$ is extended to arbitrary metric spaces, though this generalization results in an extra polylogarithmic factor in the performance guarantee.
strength_weakness: strength_weakness: The multiplicative factors in the instance optimality guarantee can be large.
questions: questions: In Theorem 2.3, the established error includes a multiplicative factor of $\log |X|$, where $|X|$ represents the size of the metric space. I would like to confirm that if $X$ is a finite $m$-dimensional space, this implies that the error scales with $m$, which can be relatively large.

I did not read the appendices.
limitations: limitations: Same as ''Questions"

Comments: title: Read the rebuttal
comment: I have read the rebuttal by the authors. I am afraid that the answers do not add much more clarity or new information to my previous assessment regarding the paper. Therefore, I'll maintain my score.