Optimal Rates for Nonparametric Density Estimation under Communication Constraints
Year: 2021

Abstract:
We consider density estimation for Besov spaces when the estimator is restricted to use only a limited number of bits about each sample. We provide a noninteractive adaptive estimator which exploits the sparsity of wavelet bases, along with a simulate-and-infer technique from parametric estimation under communication constraints. We show that our estimator is nearly rate-optimal by deriving minmax lower bounds that hold even when interactive protocols are allowed. Interestingly, while our wavelet-based estimator is almost rate-optimal for Sobolev spaces as well, it is unclear whether the standard Fourier basis, which arise naturally for those spaces, can be used to achieve the same performance.

heading: 1 Introduction
text: Estimating distributions from samples is a fundamental statistical task. Modern applications, such as those arising in federated learning or the Internet of Things (IoT), often limit access to the true data samples. One common limitation in large scale distributed systems is communication constraints, which require that each data sample must be compressed to a small number of bits.
Most prior work on communication-constrained estimation has focused on parametric problems such as Gaussian mean estimation and discrete distribution estimation. In this work, we study nonparametric density estimation under communication constraints where independent samples from an unknown distribution (whose density f lies in a function class) are distributed across players (one sample per player), and each player can only send ` bits about its sample to a central referee; the referee outputs an estimate of f based on these `-bit messages. This problem has been considered before by [5] for densities in Hölder classHsL([0, 1]) (functions supported on [0, 1] satisfying |f(x)− f(y)| ≤ L|x− y|s for every x, y ∈ [0, 1]), which is a relatively simple class for which the normalized histogram with uniform bins is known to be optimal in the centralized setting. This suggests a natural method for distributed setting – quantize each data sample into uniform bins and use the optimal estimator for distributed discrete distribution estimation. Indeed, [5] shows that this is optimal for distributed estimation of densities from the Hölder class under communication constraints. However, this simple estimator does not seem to extend to the richer Sobolev class and the most general Besov classes. In particular, the following question is largely open:
How to quantize samples to estimate densities from Besov classes under communication constraints?
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
We resolve this question both when the density belongs to a Besov class with known parameters (nonadaptive setting), and when the density belongs to a Besov class where only upper and lower bounds on parameters are known (adaptive setting). Specifically, our proposed estimators exploit the sparsity of wavelet basis for the Besov class, and use vector quantization followed by the distributed simulation technique introduced in [3] for distributed parametric estimation. We also establish lower bounds that prove the optimality of our estimators (up to logarithmic terms in the adaptive setting).

heading: 1.1 Problem setup
text: X1, . . . , Xn are independent samples from an unknown distribution with density f supported on X := [0, 1] and belonging to the Besov space B(p, q, s) (see Section 2 for details). There are n distributed users (players), with player-i having access to sample Xi. Each player can only transmit ` bits about its sample to a central server (referee) whose goal, upon observing `-bit messages from n players, is to estimate f . We consider an interactive setting, where the current player observes the messages from previous players and can use them to design its message.1 That is, in round i, player-i chooses a communication-constrained channel (randomized mapping) Wi : X → {0, 1}` as a function of prior messages Y1, . . . , Yi−1 and randomness U available to all players; it then passes Xi through Wi to generate Yi ∈ {0, 1}`. The referee observes the messages Y1, . . . , Yn and outputs an estimate f̂ of f . We term such an f̂ an (n, `)-estimate; let En,` denote the set of all (n, `)-estimates. Our goal is to design estimators that achieve the minmax expected Lr loss defined, for r ≥ 1, by
L∗r(n, `, p, q, s) := inf f̂∈En,` sup f∈B(p,q,s)
Ef [∥∥f̂ − f∥∥r
r
] . (1)
For upper bounds on L∗r(n, `, p, q, s), we consider algorithms that use the more restricted noninteractive protocols, where the channel Wi of player-i is not allowed to depend on the messages Y1, . . . , Yi−1 or on the common randomness U , but it may depend on the private randomness Ui available at player-i, where U1, . . . , Un are independent of each other and jointly of U . Noninteractive protocols are easier to implement and result in much simpler engineering for the distributed system.

heading: 1.2 Our results and techniques
text: Our first result is an information-theoretic lower bound on L∗r . Theorem 1.1. For any p, q, s, r, there exist constants C = C(p, q, s, r) > 0, α = α(p, q, s, r) > 0 such that
L∗r(n, `, p, q, s) ≥ C ·  max{n− rs 2s+1 , ( n2` )− rs2s+2 }, r ≤ (s+ 1)p, max{n− rs 2s+1 , ( n2` )− r(s−1/p+1/r)2(s−1/p)+2 (logn2`)−α}, r ∈ ((s+ 1)p, (2s+ 1)p) , max{ (
n logn
)− r(s−1/p+1/r)2(s−1/p)+1 , ( n2` )− r(s−1/p+1/r)2(s−1/p)+2 (logn2`)−α}, r ≥ (2s+ 1)p.
We emphasize that this lower bound even applies to interactive protocols as defined in Section 1.1. When the parameters p, q, s of the Besov space are known, we design a noninteractive estimator that achieves the optimal rate when r ≤ p. Theorem 1.2. For any r ≥ 1 and p, q, s with r ≤ p, there exists a constant C = C(p, q, s, r) and an (n, `)-estimate f̂ formed using a noninteractive protocol such that
sup f∈B(p,q,s)
Ef [∥∥f̂ − f∥∥r
r
] ≤ C max{n− rs 2s+1 , ( n2` )− rs2s+2 }.
We finally design an adaptive, noninteractive estimator that only requires bounds on s, and no further knowledge of p and q. Moreover, this estimator achieves (up to logarithmic factors) the optimal rate for all parameter values.
Theorem 1.3. For any N ∈ N, r ≥ 1, and p, q, s with 1/p < s < N , there exist constants C = C(p, q, s, r), α = α(p, q, s, r) and an (n, `)-estimate f̂ formed using a noninteractive protocol
1Since our lower bounds rely on general results in [2], we borrow the notation from that paper.
such that
sup f∈B(p,q,s)
Ef [∥∥f̂ − f∥∥r
r
] ≤ C logα n ·  max{n− rs 2s+1 , ( n2` )− rs2s+2 }, r ≤ (s+ 1)p, max{n− rs 2s+1 , ( n2` )− r(s−1/p+1/r)2(s−1/p)+2 }, r ∈ ((s+ 1)p, (2s+ 1)p) ,
max{n− r(s−1/p+1/r) 2(s−1/p)+1 , ( n2` )− r(s−1/p+1/r)2(s−1/p)+2 }, r ≥ (2s+ 1)p,
where the protocol only requires knowledge of N (an upper bound on s).
In summary, for all r ≥ 1, the minmax Lr loss of estimating B(p, q, s) (up to logarithmic factors) is
L∗r(n, `, p, q, s)  max{n− rs 2s+1 , ( n2` )− rs2s+2 } r ≤ (s+ 1)p, max{n− rs 2s+1 , ( n2` )− r(s−1/p+1/r)2(s−1/p)+2 } r ∈ ((s+ 1)p, (2s+ 1)p) ,
max{n− r(s−1/p+1/r) 2(s−1/p)+1 , ( n2` )− r(s−1/p+1/r)2(s−1/p)+2 } r ≥ (2s+ 1)p. (2)
It it worth noting that the first term of the maximum in all cases is the standard, unconstrained nonparametric rate (cf. [10], or the discussion below), while the second term reflects the convergence slowdown due to the communication constraints. The effect of communication constraints disappears when ` is sufficiently large. In particular, we get back the centralized rates when ` satisfies
` ≥  ( 1 2s+1 ) logn, if r ≤ (s+ 1)p,( 2s(1−1/r) (2s+1)(s−1/p+1/r) − 1 2s+1 ) logn, if r ∈ ((s+ 1)p, (2s+ 1)p) ,(
1 2(s−1/p)+1
) logn, if r ≥ (2s+ 1)p.
For the standard case of L2 loss, with, say p ≥ 2, the minmax rate becomes the more interpretable quantity
L∗2(n, `, p, q, s) max{n− 2s 2s+1 , ( n2` )− 2s2s+2 },
where we see that the `-bit communication constraint reduces the exponent of the convergence rate from 2s2s+1 to 2s 2s+2 . Note that the difference becomes less perceptible as ` grows, or s tends to∞. Finally, from (2) we observe qualitative changes at r = (s + 1)p and r = (2s + 1)p, where the rate exponent changes slope. This phenomenon, whose analogue is observed in the unconstrained setting [10] as well as under local privacy constraints [7], is sometimes referred to as an elbow effect.
Quantize, simulate and infer. A conceptually simple technique for distributed estimation under communication constraints (“simulate-and-infer”) was proposed in [3], which uses communication to simulate samples from the unknown distribution, and provides an optimal rate estimator for discrete distribution estimation under communication constraints. A natural approach for nonparametric estimation would be to quantize the samples to the available number of bits (`) and use this quantized sample to estimate the distribution. However, it is unclear if this approach gives optimal rates. Instead, in our approach, we quantize without inducing any “loss of information.” Specifically, we form an approximately sufficient statistic (based on wavelets) that can be represented using finite bits and does not result in rate-loss. The number of bits could still be more than `, and therefore, we use simulate-and-infer to generate samples from the statistic. Thus, loss of information due to communication constraints only happens in the last step, when we use multiple samples to simulate a sample from the sufficient statistic; the quantization part is just for efficient finite representation.
Sobolev spaces and Fourier bases. A first approach that we tried for Sobolev spaces was to use the Fourier basis, a natural choice for a Sobolev space. However, all our attempts led to a suboptimal performance either in the dependence on ` (we were not able to get the exponential 2−` dependence) or the exponent of n (i.e., when we tried to get 2−` dependence, this resulted in a suboptimal exponent of n). Somewhat surprisingly, the more general wavelet-based approach described above gives us tight bounds for Sobolev space as well, since Sobolev space S(β) = B(2, 2, β). Thus, it seems that it is necessary to use the wavelet representation even for Sobolev spaces to get an appropriately “small” statistic – sparsity of wavelets is very useful for inference under communication constraints.
Organization. After discussing prior works (Section 1.3) and preliminaries (Section 2), we describe our estimation algorithms in Section 3. In Section 3.1, we present the nonadaptive single-level estimator (comprising Algorithm 2 and Algorithm 3), which achieves guarantees of Theorem 1.2
(proof in Appendix B). In Section 3.2, we present the adaptive multi-level estimator (comprising Algorithm 4 and Algorithm 5), which achieves guarantees of Theorem 1.3 (proof in Appendix C). Finally, in Section 4, we briefly discuss our approach to proving lower bounds of Theorem 1.1, with details deferred to Appendix D. Appendices can be found in the Supplementary.

heading: 1.3 Prior work
text: The seminal work of Donoho, Johnstone, Kerkyacharian, and Picard [10] proposed wavelet estimators for Besov class in the centralized setting, and showed that these estimators achieve near-optimal rates of convergence (up to logarithmic factors),
L∗r(n,∞, p, q, s) { n− rs 2s+1 , r < (2s+ 1)p,
n− r(s−1/p+1/r) 2(s−1/p)+1 , r ≥ (2s+ 1)p. (3)
Their results highlight the fact that linear estimators are inherently suboptimal for estimation with respect to Lr losses, when r is large; that is, some nonlinearity in the estimator is required to achieve optimal rates for r > p. In particular, they show that nonlinearity in the form of thresholding achieves optimal rates for r > p (see Section 2 below for details). Further, they use thresholding to design adaptive estimators that achieve near-optimal rates. These minmax rates exhibit the aforementioned elbow effect, where the error exponent is only piecewise linear, and changes slope at r = (2s+ 1)p. We refer the reader to [10] for a further discussion of these phenomena.
Butucea, Dubois, Kroll, and Saumard [7] recently extended these ideas to obtain near-rate optimal estimators for Besov spaces under local differential privacy constraints. Their adaptive estimator, as well as the information-theoretic lower bounds they establish, show that similar phenomena occur in the context of locally private nonparametric estimation. Our work, specifically the analysis of our adaptive estimator, draws upon some of the ideas of [7], with some crucial differences. In particular, the key ideas underlying our estimators – the wavelet-induced sparsity (Claim 3.1), the use of distributed simulation, and vector quantization – are neither present in nor applicable to the setting of [7] (where the introduction of random noise to ensure differential privacy effectively removes wavelet sparsity). Furthermore, our lower bounds even apply to interactive protocols, unlike the lower bounds from [7] which are restricted to the noninteractive setting.
In summary, our paper is the first to derive the counterpart of the nonparametric estimation results of [10, 7] under communication constraints, and shows that the analogue of the phenomena observed in [10] holds in the communication-constrained setting.
Other work on distributed estimation. We briefly discuss the related literature on distributed (communication-constrained) estimation problems. [18, 15, 16, 8] have studied the problem of distributed nonparametric function estimation (regression) under a Gaussian white noise model in a noninteractive setting with n players, where each player observes an independent copy of the stochastic process dY (t) = f(t)dt+ (1/ √ N)dW (t), 0 ≤ t ≤ 1. Here W (t) is the standard Wiener process, and f is the function to be estimated. [18] derive minimax rates for f in Sobolev space under L2 loss, where each player can send at most ` bits. [15] derive minimax rates for f in the Besov space B(2,∞, s) (“Sobolev type”) under L2 loss, and f in B(∞,∞, s) (“Hölder type”) under L∞ loss, where each player can send at most ` (assumed to be at least logN ) bits on average. Further, the paper proposes near-optimal adaptive estimators (based on Lepski’s method) that adapt to the smoothness parameter s, provided that s ∈ [smin, smax), where smin depends on n, ` and smax can be arbitrary. [16] further study the problem of adaptivity for B(2,∞, s) under L2 loss and B(∞,∞, s) under L∞ loss, and answer the question of whether it is possible to design adaptive estimators (adapting to the smoothness parameter s) that attain centralized minimax rates while also having the expected communication budget nearly the same as that of a minimax optimal distributed estimator that knows s. The paper shows that this is possible for B(2,∞, s) under L2 loss provided that s is below a certain threshold, and is impossible for B(∞,∞, s) under L∞ loss. [8] derive minimax rates for f in Besov spaces B(p, q, s) (for p ≥ 2) under L2 loss, where each player can send at most ` bits on average. In addition, they study the problem of adaptivity and characterize the minimax communication budget of adaptive estimators (adapting to parameters p ≥ 2, q > 1, s > 0) that achieve centralized rates. The adaptive estimator proposed in this paper is based on thresholding, where the thresholding is done locally by each player.
Interestingly, when the communication budget ` is insufficient to achieve centralized rates, the minimax rates in these distributed nonparametric function estimation problems decay polynomially in
`, which is in contrast to the minimax rates we obtain for distributed nonparametric density estimation problem, where the decay is exponential in `.
For higher dimensions, [13] consider two-party estimation of Hölder smooth density functions in d dimensions, where the two parties each observe a subset of coordinates and must communicate to estimate the density. They show that interactivity strictly helps over one-way communication.
We now discuss related works on distributed parametric estimation problems. [5] establish lower bounds on parametric density estimation, and on some restricted nonparametric families (Hölder classes) by bounding the Fisher information. [3] obtain upper and lower bounds for discrete distribution estimation; our algorithms leverage the concept of distributed simulation (“simulate-and-infer”) introduced in that context. [12] derive lower bounds for various parametric estimation tasks, including discrete distributions and continuous (parametric) families such as high-dimensional Gaussians (including the sparse case). [2], building on [1] (which focused on learning and testing discrete distributions), developed a general technique to prove estimation lower bounds for parametric families; our lower bounds rely on their framework, by suitably extending it to handle the nonparametric case.
We note that there are other approaches for establishing lower bounds under communication constraints such as the early works [17], [14]; [11, 6], where bounds for specific inference problems under communication constraints were obtained; and [4], where Cramér–Rao bounds for this setting were developed. We found the general approach of [2] best fits our specific application, where we needed to handle interactive communication as well as a nonuniform prior on the parameter in the lower bound construction.

heading: 2 Preliminaries
text: Given two integersm ≤ n, we write Jm,nK for the set {m,m+1, . . . , n} and JnK for J1, nK. For two sequences or functions (an)n, (bn)n, we write an . bn if there exists a constant C > 0 (independent of n) such that an ≤ Cbn for all n, and an bn if both an . bn and an & bn. For a function g, supp(g) denotes the support of g. Let B(p, q, s) be the Besov space with parameters p, q, s, where 1 ≤ p, q ≤ ∞ and s > 0. We will write φ, ψ ∈ L2(R), respectively, for the father and mother wavelets generating the basis of the Besov space, and ‖f‖pqs for the Besov norm of f ∈ B(p, q, s). We refer the reader to Appendix A in Supplementary for details on wavelets, Besov spaces and Besov norm.
Assumptions on the density and wavelets. We make the following assumptions on the density f :
1. f is compactly supported: without loss of generality, supp(f) ⊆ [0, 1]. 2. Besov norm of f is bounded: without loss of generality, ‖f‖pqs ≤ 1.
Our algorithm works with any father and mother wavelets φ and ψ satisfying the following conditions:
1. φ and ψ are N -regular, where N > s, and 2. supp(φ), supp(ψ) ⊆ [−A,A] for some integer A > 0 (which may depend on N ).
As a concrete example, Daubechies’ family of wavelets [9] satisfies these assumptions.
Density estimation in centralized setting. In the centralized setting, X1, . . . , Xn from an unknown density f ∈ B(p, q, s) are accessible to the estimator. Let the wavelet expansion of f be (see Appendix A.1 for details)
f = ∑ k∈Z α0,kφ0,k + ∑ j≥0 ∑ k∈Z βj,kψj,k (4)
where φj,k(x) = 2j/2φ(2jx − k), ψj,k(x) = 2j/2ψ(2jx − k). The wavelet basis satisfies the property that, for any L,H ∈ Z with H ≥ L, we have ∑ k∈Z αL,kφL,k + ∑H−1 j=L ∑ k∈Z βj,kψj,k =∑
k∈Z αH,kφH,k (Fact A.1 in Supplementary). Note that for a given j, k, α̂j,k := 1 n ∑n i=1 φj,k(Xi)
is an unbiased estimate of αj,k. Thus, for some H ∈ Z+, an estimate of f is
f̂lin = ∑ k∈Z α̂H,kφH,k, α̂H,k := 1 n n∑ i=1 φH,k(Xi), (5)
where H is chosen depending on n and parameters p, q, s to minimize the (worst-case) Lr loss. This simple estimator (with appropriate choice of H) is rate-optimal when 1 ≤ r ≤ p, but is sub-optimal when r > p [10]. Moreover, setting H requires knowing the Besov parameters p, q, s, which renders this estimator nonadaptive. The main contribution of [10] was to demonstrate that thresholding leads to estimators that are (i) near-optimal for every r ≥ 1; (ii) adaptive, in the sense that the estimator does not use the values of parameters p, q, s as long as 1/p < s < N for some N ∈ N. For a given L,H ∈ Z+, L ≤ H , a thresholded estimator outputs the estimate
f̂thresh = ∑ k∈Z α̂L,kφL,k + H∑ j=L ∑ k∈Z β̃j,kψj,k, β̃j,k = β̂j,k1{|β̂j,k|≥tj} (6)
where α̂L,k := 1n ∑n i=1 φL,k(Xi), β̂j,k = 1 n ∑n i=1 ψj,k(Xi), and tj is a fixed threshold proportional
to √ j/n; here, L,H depend on n, but not on parameters p, q, s. Our proposed estimators draws upon these classical estimators.

heading: 3 Estimation Algorithms
text: We propose algorithms for density estimation under communication constraints that achieve optimal/near-optimal performance in terms of n (number of players) and ` (number of bits each player can send). Designing a density estimator in the communication-constrained setting consists of: (i) specifying the sample-dependent `-bit message that a player sends to the referee; (ii) specifying the density estimate that the referee outputs based on the `-bit messages from the n players. As in the unconstrained setting, we estimate f by estimating its wavelet coefficients.
Our estimators consist of three ingredients: wavelet-induced sparsity, vector quantization, and distributed simulation.
(i) Wavelet-induced sparsity. Let the wavelet expansion of the density function f be given by (4). For a given J ∈ Z+, partition the interval [0, 1] into 2J uniform bins as
[0, 1] = 2J−1⋃ t=0 E (J) t where E (J) t := {[ t2−J , (t+ 1)2−J ) if t ∈ J0, 2J − 2K,[ 1− 2−J , 1 ] if t = 2J − 1. (7)
For a bin E(J)t , t ∈ J0, 2J − 1K, let A(J)t := { k ∈ Z : E(J)t ∩ supp(φJ,k) is non-empty } ; (8)
B(J)t := { k ∈ Z : E(J)t ∩ supp(ψJ,k) is non-empty } . (9)
That is, for x ∈ E(J)t , we have φJ,k(x) = 0 for k /∈ A (J) t , and ψJ,k(x) = 0 for k /∈ B (J) t . By “wavelet-induced sparsity,” we mean the following:
Claim 3.1. Let [0, 1] = ⋃2J−1 t=0 E (J) t as in (7). Then, for each t ∈ J0, 2J − 1K,
|A(J)t | ≤ 2(A+ 2), |B (J) t | ≤ 2(A+ 2),
where A is the assumed bound for points in the support of φ and ψ.
The claim follows from the observation that φJ,k (resp., ψJ,k) is obtained by translating φJ,0 (resp., ψJ,0) in steps of size 2−J , and that supp(φJ,k) ⊆ [−A2−j , A2−j ] (resp., supp(ψJ,k) ⊆ [−A2−j , A2−j ]). (ii) Vector quantization. Consider the problem of designing a randomized algorithm that takes as input an arbitrary x ∈ Rd satisfying ‖x‖∞ ≤ B, and outputs a random vector Q(x) ∈ Rd chosen from an alphabet of finite cardinality, such that E[Q(x)] = x. Our vector quantization algorithm (Algorithm 1) achieves this, and is based on the following idea: Let P be a convex polytope with vertices {v1, v2, . . .} such that { x ∈ Rd : ‖x‖∞ ≤ B } ⊆ P . Given x (with ‖x‖∞ ≤ B), express x
as convex combination of vertices of P (say, x = ∑ i θivi) and output a random vertex V , where V = vi with probability θi. Clearly, E[V ] = x.
Specifically, Algorithm 1 uses the polytope P = PV formed by the vertex set V = {±(Bd)e1, . . . ,±(Bd)ed}, where ei is the i-th standard basis vector (i.e., P is the `1 ball of radius Bd). Note that |V| = 2d and that { x ∈ Rd : ‖x‖∞ ≤ B } ⊆ PV . This leads to the following claim.
Claim 3.2. Given x ∈ Rd with ‖x‖∞ ≤ B as input, Algorithm 1 outputs a random variable Q(x) ∈ V that is an unbiased estimate of x, with |V| = 2d.
Remark. A more direct approach to quantization would be to do it coordinate-wise, i.e., quantize (independently) each coordinate to {−B,B} with appropriate probability to make it unbiased. This can equivalently be seen as quantizing the vector using the `∞ ball (of radius B) as the polytope. Here, the alphabet size becomes 2d instead of 2d in Algorithm 1; but, on the plus side, the coordinatewise variance of the quantized vector becomes ≈ B2, instead of ≈ (Bd)2 in Algorithm 1. In our estimators, we will be quantizing vectors of constant length (d), so these dependencies on d do not affect the rate (up to constants).
Algorithm 1 Vector quantization Let V = {±(Bd)e1, . . . ,±(Bd)ed}. Label the vectors in V as v1, . . . , v2d. Input: x ∈ Rd with ‖x‖∞ ≤ B.
1: Write x as convex combination of vectors in V: x = ∑2d i=1 θivi. 2: Choose I ∈ {1, . . . , 2d} randomly where I = i with probability θi and return Q(x) = vI .
(iii) Distributed simulation. The problem of distributed simulation is the following: There are n players, each having an i.i.d. sample from an unknown d-ary distribution p. Each player can only send ` bits to a central referee, where ` < log d. Can the referee simulate i.i.d. samples from p using `-bit messages from the players? [3] proposed a noninteractive communication protocol, using which the referee can simulate one sample from p using `-bit messages from O(d/2`) players. Moreover, the protocol is deterministic at the players, and only requires private randomness at the referee. Theorem 3.3 ([3]). For any ` ≥ 1, the simulation protocol of [3], denoted DISTRSIM`, lets the referee simulate Ω(n2`/d) i.i.d. samples from an unknown d-ary probability distribution p using `-bit messages from n players, where each player holds an independent sample from p.
Combining ideas. We now discuss how the three ideas come together. To mimic the classical estimator (5), a player with sample X would ideally like to communicate {φH,k(X)}k∈Z, but cannot do so due to communication constraints. Wavelet-induced sparsity (Claim 3.1) ensures that communicating the bin (out of 2H possible bins) in which X lies is tantamount to identifying the set of at most d := 2(A+2) indices k for which φH,k(X) is possibly non-zero. Moreover, the player can quantize (unbiasedly) the vector containing values of φH,k(X) at these indices using Algorithm 1, whose output is one of 2d possibilities (Claim 3.2). Thus, overall, using an alphabet of size at most 2H × (2d) = O(2H), a player can communicate an unbiased estimate of {φH,k(X)}k∈Z. It can be shown that a density estimate based on these unbiased estimates from n players still achieves centralized minmax rates (up to constants). However, if 2` < 4(A+ 2)2H , the players cannot send these estimates directly to the referee. In this case, the players and the referee use the distributed simulation protocol DISTRSIM` (Theorem 3.3), which, effectively, enables the referee to simulate O(n2`/2H) i.i.d. realizations of unbiased estimates of {φH,k(X)}k∈Z. The referee can now output a density estimate based on these simulated estimates. The degradation in minmax rates under communication constraints is due to the fact that the referee has only O(n2`/2H) realizations of unbiased estimates of {φH,k(X)}k, instead of n. We now give details of the idea outlined above. The resulting estimator (single-level estimator) is a communication-constrained version of the classical estimator given in (5). We then describe an adaptive estimator (multi-level estimator), which is a communication-constrained version of the classical adaptive estimator (6).

heading: 3.1 Single-level estimator
text: The n players and the referee agree beforehand on the following: wavelet functions φ, ψ; H ∈ Z+; partition [0, 1] = ⋃2H−1 t=0 E (H) t as in (7); collections of indices A (H) t , t ∈ J0, 2H − 1K as in (8). For every t, the indices in A(H)t are arranged in ascending order.
Player’s side. Each player carries out two broad steps: (i) quantization; (ii) simulation.
Algorithm 2 Single-level estimator (Players) Input: Player-i has input Xi, i ∈ JnK.
1: for i = 1, . . . , n do 2: Player-i computes Zi := (Bi, Q(Vi)), where: (i) Bi is the bin in which Xi lies; (ii) Q(Vi) is an unbiased quantization of the vector Vi := { 2−H/2φH,k(Xi) } k∈A(H)
Bi . . Quantization 3: Player-i computes `-bit message Yi corresponding to Zi as per DISTRSIM` (Theorem 3.3),
and sends it to the referee. . Simulation
The scaling by 2−H/2 in the definition of Vi (line 2 in Algorithm 2) ensures that ‖Vi‖∞ ≤ ‖φ‖∞, which is a constant. This enables the use of Algorithm 1 to compute quantization of Vi. Overall, computing Zi = (Bi, Q(Vi)) involves two quantizations: Bi can be seen as a quantized version of Xi ∈ [0, 1]; Q(Vi) is a quantized version of {φH,k(Xi)}k. Moreover, for each i ∈ JnK, Zi ∈ Z(H), where Z(H) := J0, 2H − 1K × {±Be1, . . . ,±Bed} (with d ≤ 2(A + 2), by Claim 3.1), so that∣∣Z(H)∣∣ ≤ 4(A+ 2) 2H = O(2H). Thus, Z1, . . . , Zn are i.i.d. samples (since X1, . . . , Xn are i.i.d.) from a
∣∣Z(H)∣∣-ary distribution (call it pZ(H) ) distributed across n players, where
∣∣Z(H)∣∣ = O(2H). Since a player can send only ` bits, player-i cannot send Zi directly if 2` <
∣∣Z(H)∣∣. In this case, player-i computes an `-bit message Yi according to the distributed simulation protocol DISTRSIM`, and sends Yi to the referee.
Referee’s side. The referee, using the simulated i.i.d. samples from pZ(H) , computes the density estimate similar to the classical estimate (5). This is possible because the m = O(n2`/2H) simulated samples are, essentially, i.i.d. realizations of unbiased quantization of {φH,k(X)}k∈Z.
Algorithm 3 Single-level estimator (Referee) Input: Y1, . . . , Yn (`-bit messages from n players).
1: From Y1, . . . , Yn, referee obtains m = O(n2`/ ∣∣Z(H)∣∣) = O(n2`/2H) i.i.d. samples
Z ′1, . . . , Z ′ m ∼ pZ(H) as per DISTRSIM`, where Z ′i = (B′i, Q′i) ∈ Z(H).
2: for i = 1, . . . ,m do 3: Referee computes
φ̂ (i) H,k :=
{ 2H/2 Q′i(k) if k ∈ A (H) B′
i 0 otherwise, (10)
where Q′i(k) is the entry in Q′i corresponding to index k ∈ A (H) B′ i . . Scaling by 2H/2 is to
negate the scaling by 2−H/2 used in definition of Vi on the players’ side. 4: Referee outputs density estimate
f̂ = ∑ k∈Z α̂H,kφH,k, where α̂H,k = 1 m m∑ i=1 φ̂ (i) H,k, k ∈ Z. (11)
Result. For H such that 2H min{(n2`) 1 2s+2 , n 1
2s+1 }, the single-level estimator recovers the guarantees in Theorem 1.2 (see Appendix B in Supplementary). The estimator is nonadaptive because setting H requires knowing Besov parameter s. Further, note that the estimator is indeed noninteractive, as player-i’s message Yi does not depend on messages Y1, . . . , Yi−1.

heading: 3.2 Multi-level estimator: An adaptive density estimator
text: The key observation in designing our multi-level estimator is that different coefficients need to be recovered with different accuracy. We enable this by dividing players into groups for estimating different coefficients, and using a different level of quantization for each group. This is in contrast to simply mimicking the classical adaptive estimator (6), which would suggest that a player with sample X should quantize and communicate information about {φL,k(X)}k , { {ψJ,k(X)}k } J∈JL,HK. Instead, we do the following: Divide n players into H − L + 1 groups of equal size (so, each
group has n′ = nH−L+1 players). Label the groups L,L+ 1, . . . ,H . Players in group-L only focus on {φL,k(X)}k , {ψL,k(X)}k. Players in group-J , J ∈ JL + 1, HK, only focus on {ψJ,k(X)}k. Moreover, players in group-J , J ∈ JL,HK quantize their sample X using 2J uniform bins. As before, by Claim 3.1, this is tantamount to identifying at most a constant number of indices for which the wavelet function evaluates to a non-zero value (since players in group-J only consider φJ,k or ψJ,k). The player then quantizes the vector containing these values using Algorithm 1, before using distributed simulation.
The n players and the referee agree beforehand on the following: wavelet functions φ, ψ; L,H ∈ Z+; division of players into H − L+ 1 groups. Further, for each J ∈ JL,HK, the n′ players in group-J and the referee agree on the following: partition [0, 1] = ⋃2J−1 t=0 E (J) t as in (7); collection of indices A(J)t ,B (J) t , t ∈ J0, 2J − 1K, as in (8), (9). For every J, t, the indices in A (J) t ,B (J) t are arranged in ascending order.
Player’s side. Label players in group-J as (1, J), . . . , (n′, J). We denote by Xi,J the sample with player-(i, J). Essentially, players in group-J run quantization and simulation steps as in the single-level algorithm (Algorithm 2), with H replaced by J .
Algorithm 4 Multi-level estimator (Players) Input: Player-(i, J) has input Xi,J , i ∈ Jn′K, J ∈ JL,HK (where n′ = nH−L+1 ).
1: for J = L,L+ 1, . . . ,H do 2: for i = 1, . . . , n′ do 3: Player-(i, J) computes Zi,J = (Bi,J , Q(Vi,J)), where: (i)Bi,J is the bin (out of 2J bins)
in which Xi,J lies; (ii) Q(Vi,J) is an unbiased quantization of the vector Vi,J , where
Vi,J :=  { 2−J/2ψJ,k(Xi,J) } k∈B(J) Bi,J if J ∈ JL+ 1, HK,{ 2−L/2φL,k(Xi,L) } k∈A(L)
Bi,L
⊕ { 2−L/2ψL,k(Xi,L) } k∈B(L)
Bi,L
if J = L.
(⊕ denotes concatenation of two vectors.) . Quantization 4: Player-(i, J) computes `-bit message Yi,J corresponding to Zi,J as per DISTRSIM`, and
sends it to the referee. . Simulation Referee’s side. For players in group-J , Zi,J ∈ Z(J), where ∣∣Z(J)∣∣ = O(2J). Thus, after distributed simulation, referee obtains mJ = O(n′2`/2J) = O(n2`/(H − L + 1)2J) samples from players of group-J . Note that, higher the J , fewer the simulated samples; this dependence on J of the number of samples available with the referee is one of the major differences between the classical and the distributed setting. Finally, using the simulated samples from players of every group, referee computes a density estimate similar to the adaptive classical estimator (6), with threshold value tJ = κ √ J/mJ , for a constant κ.
Result. For L,H satisfying 2L min{(n2`) 1 2(N+1)+2 , n 1 2(N+1)+1 } and 2H min{ √ n2`
log(n2`) , n logn}, the multi-level estimator yields the guarantees in Theorem 1.3 (see Appendix C in Supplementary) as long as s ∈ (1/p,N + 1) (recall that N is the regularity of the wavelet basis). Since L,H do not depend on specific Besov parameters, the estimator is adaptive. Moreover, it is noninteractive.

heading: 4 Lower Bounds
text: We conclude with a description of our information-theoretic lower bounds (Theorem 1.1) for the minimax loss L∗r(n, `, p, q, s), which applies to the broader class of interactive protocols (recall that our matching upper bounds are obtained by noninteractive ones); the details can be found in Appendix D in Supplementary. To derive lower bounds, we consider a family of probability distributions P parameterized by {−1, 1}d for some d ∈ Z+; that is, P = {pz : z ∈ {−1, 1}d}, where pz has density fz . Moreover, we specify a prior π on Z = (Z1, . . . , Zd) ∈ {−1, 1}d, defined as Zi ∼ Rademacher(τ) independently for each i ∈ [d], for some τ ∈ (0, 1/2]. We then consider the following scenario:
Algorithm 5 Multi-level algorithm (Referee) Input: {Yi,J}i∈Jn′K,J∈JL,HK (`-bit messages from n players).
1: for J = L,L+ 1, . . . ,H do 2: From Y1,J , . . . , Yn′,J , referee obtains mJ = O(n2`/(H − L + 1)2J) i.i.d. samples Z ′1,J , . . . , Z ′ mJ ,J
∼ pZ(J) as per DISTRSIM`, where Z ′i,J = (B′i,J , Q′i,J) ∈ Z(J). 3: for i = 1, . . . ,mJ do 4: if J = L then
5: Referee computes { φ̂
(i) L,k } k∈Z as φ̂(i)L,k := { 2L/2 Q′i,L(k) if k ∈ A (L) B′ i,L
0 otherwise. 6: Referee computes { ψ̂
(i) J,k } k∈Z as ψ̂(i)J,k := { 2J/2 Q′i,J(k) if k ∈ B (J) B′ i,J
0 otherwise. 7: Referee outputs density estimate
f̂ = ∑ k α̂L,kφL,k + H∑ J=L ∑ k β̃J,kψJ,k, (12)
where α̂L,k = 1mL ∑mL i=1 φ̂ (i) L,k, β̂J,k = 1 mJ ∑mJ i=1 ψ̂ (i) J,k, β̃J,k = β̂J,k1{|β̂J,k|≥tJ :=κ√J/mJ}.
For Z ∼ π, let X1, . . . , Xn be i.i.d. samples from pZ distributed across n players. Let Y1, . . . , Yn be `-bit messages sent by the players (possibly interactively) to the referee. Denote by pY n+i (resp. pY n
−i ) the joint distribution of Y1, . . . , Yn, given Zi = 1 (resp. Zi = −1). That is,
pY n +i = 1 τ ∑ z:zi=1 π(z)pY n z , pY n −i = 1 1− τ ∑ z:zi=−1 π(z)pY n z , (13)
where pY nz is the joint distribution of Y1, . . . , Yn, given Z = z.
In this scenario, we analyze the “average discrepancy” 1d ∑d i=1 dTV ( pY n−i ,pY n +i ) , where dTV(p,q) denotes the total variation distance between p and q. On the one hand, a result from [2] gives us an upper bound on this average discrepancy as a function of n and ` which holds for any interactive protocol generating Y1, . . . , Yn (Theorem D.3 in Supplementary). On the other hand, we derive a lower bound on average discrepancy (as a function of the error rate ε) as follows: Consider a communication-constrained density estimation algorithm (possibly interactive) which outputs f̂ satisfying supf∈B(p,q,s) Ef [ ‖f̂ − f‖rr ] ≤ εr. We show that one can use the messages Y1, . . . , Yn generated by this algorithm to solve, for each i ∈ JdK, the binary hypothesis testing problem of deciding whether Zi = 1 or Zi = −1. This, in turn, implies a lower bound on 1 d ∑d i=1 dTV ( pY n−i ,pY n +i ) . Putting together the upper and lower bounds on average discrepancy gives us a lower bound on ε.
The parameterized family of distributions P we consider is constructed as follows: Let f0 be a function supported on [0, 1]. Let I1, . . . , Id ⊆ [0, 1] be mutually disjoint intervals of equal length. Let ψi be a “bump” function supported on interval Ii, where ψi’s are all translations of the same bump function. Then, for z = (z1, . . . , zd) ∈ {−1, 1}d, we define pz to be a probability distribution with density fz , defined as the “baseline” f0 perturbed by adding (a rescaling of) the bump ψi according to the value of zi. In more detail, to get the desired lower bounds, we distinguish two cases depending on whether r < (s+ 1)p, and construct two families of distributions: P1 (when r < (s+ 1)p) and P2 (when r ≥ (s+ 1)p). For P1, we use a uniform prior on Z = (Z1, . . . , Zd), i.e., Z has independent Rademacher(1/2) coordinates, and set fz = f0 + γ ∑d i=1 ziψi for some suitably small parameter γ > 0. That is, the baseline density f0 has disjoint bumps, which are either ψi or −ψi depending on the value of zi. (See Appendix D.3 for details.) For P2, we use a non-uniform (“sparse”) prior on Z, where Z has independent Rademacher(1/d) coordinates, and set fz = f0 + γ ∑d i=1(1 + zi)ψi (so that bump ψi only appears if Zi = 1). (See Appendix D.4 for details.) Applying to these constructions the method described above allows us to derive the lower bounds of Theorem 1.1.
Acknowledgments and Disclosure of Funding
Aditya Vikram Singh is supported by Prime Minister’s Research Fellowship. Jayadev Acharya is supported by NSF-CCF- 1846300 (CAREER), NSF-CCF-1815893, and a Google Faculty Fellowship.

Review summaries:
This problem of estimating a density, assumed to lie in a Besov class, under communication constraints in which each sample must be encoded using a fixed number $\ell$ of bits before computing the density estimate. Both the non-adaptive and adaptive cases are considered, with matching upper and lower bounds (up to polylogarithmic factors) obtained in both cases. The matching bounds also indicate that the rate is the same for interactive and non-interactive protocols. The proposed methods utilizes 3 main ingredients: 1. The fact that at most a constant number of wavelets at resolution J are non-zero on any interval of width $2^{-J}$. 2. A particular randomized scheme that quantizes the values of the non-zero wavelets at a sample point without introducing bias. 3. The DISTRSIM protocol proposed by [3], which maps a quantized sample to a smaller IID sample of higher resolution.

The paper proposed a method for distributed density estimation with the unknown density from a Besov space. The estimation is based on local observations of independent samples drawn from a common distribution. The agents with local observations are connected via bit-constrained channels to a central unit ("referee"). The proposed method is based on a novel quantization strategy that exploits the properties of a wavelet expansion of the unknown density. The method is implemented as two algorithms, one which requires knowledge of Besov space parameter and another one which does not require precise Besov parameters.  

The paper addresses learning densities in Bezov space under communication constraints, primarily in a non-interactive setting. In fact, Theorems 1.1 and 1.3 together suggest that interaction may not help much in achieving a better performance. The proposed algorithms leverage wavelet-induced sparsity, vector quantization, and distributed (sample) simulation, but not the standard Fourier basis.