Density estimation for shift-invariant multidimensional distributions
Anindya De∗
Northwestern University
anindya@eecs.northwestern.edu

Philip M. Long
Googleplong@google.com

Rocco A. Servedio†
Columbia University
rocco@cs.columbia.edu

Abstract
We study density estimation for classes of shift-invariant distributions over Rd . A multidimensional
distribution is “shift-invariant” if, roughly speaking, it is close in total variation distance to a small shift of
it in any direction. Shift-invariance relaxes smoothness assumptions commonly used in non-parametricdensity estimation to allow jump discontinuities. The different classes of distributions that we consider
correspond to different rates of tail decay.
For each such class we give an efficient algorithm that learns any distribution in the class fromindependent samples with respect to total variation distance. As a special case of our general result,
we show that d-dimensional shift-invariant distributions which satisfy an exponential tail bound can belearned to total variation distance error ε using Õd (1/εd+2 ) examples and Õd (1/ε2d+2 ) time. Thisimplies that, for constant d, multivariate log-concave distributions can be learned in Õd (1/ε2d+2 ) timeusing Õd (1/εd+2 ) samples, answering a question of [DKS16b]. All of our results extend to a model ofnoise-tolerant density estimation using Huber’s contamination model, in which the target distribution to
be learned is a (1 − ε, ε) mixture of some unknown distribution in the class with some other arbitrary
and unknown distribution, and the learning algorithm must output a hypothesis distribution with total
variation distance error O(ε) from the target distribution. We show that our general results are close tobest possible by proving a simple Ω 1/εd information-theoretic lower bound on sample complexity
even for learning bounded distributions that are shift-invariant.

1

Introduction

In multidimensional density estimation, an algorithm has access to independent draws from an unknown
target probability distribution over Rd , which is typically assumed to belong to or be close to some classof “nice” distributions. The goal is to output a hypothesis distribution which with high probability is close
to the target distribution. A number of different distance measures can be used to capture the notion ofcloseness; in this work we use the total variation distance (also known as the “statistical distance” andequivalent to the L1 distance). This is a well studied framework which has been investigated in detail, seee.g. the books [DG85, DL12].
∗†

Supported by NSF grant CCF-1814706
Supported by NSF grants CCF-1319788 and CCF-1420349

1

Multidimensional density estimation is typically attacked in one of two ways. In the first general approach a parameterized hypothesis class is chosen, and a setting of parameters is chosen based on the observed data points. This approach is justified given the belief that the parameterized class contains a good
approximation to the distribution generating the data, or even that the parameterized class actually contains
the target distribution. See [Das99, KMV10, MV10] for some well-known multidimensional distribution
learning results in this line.
In the second general approach a hypothesis distribution is constructed by “smoothing” the empirical
distribution with a kernel function. This approach is justified by the belief that the target distribution satisfies
some smoothness assumptions, and is more appropriate when studying distributions that do not have aparametric representation. The current paper falls within this second strand.
The most popular smoothness assumption is that the distribution has a density that belongs to a Sobolev
space [Sob63, BC91, HK92, DL12]. The simplest Sobolev space used in this context corresponds to having abound on the average of the partial first “weak derivatives” of the density; other Sobolev spaces correspondto bounding additional derivatives. A drawback of this approach is that it does not apply to distributions
whose densities have jump discontinuities. Such jump discontinuities can arise in various applications, forexample, when objects under analysis must satisfy hard constraints.
To address this, some authors have used the weaker assumption that the density belongs to a Besov
space [Bes59, DS93, Mas97, WN07, ADLS17]. In the simplest case, this allows jump discontinuities aslong as the function does not change very fast on average. The precise definition, which is quite technical(see [DS93]), makes reference to the effect on a distribution of shifting the domain by a small amount.
The densities we consider. In this paper we analyze a clean and simple smoothness assumption, which is acontinuous analog of the notion of shift-invariance that has recently been used for analyzing the learnabilityof various types of discrete distributions [BX99, DDO+ 13, DLS18]. The assumption is based on the shiftinvariance of f in direction v at scale κ, which, for a density f over Rd , a unit vector v ∈ Rd , and a positive
real value κ, we define to beZ
def 1SI(f, v, κ) = · supf (x + κ0 v) − f (x) dx.
κ κ0 ∈[0,κ] Rd
We define the quantity SI(f, κ) to be the worst case of SI(f, v, κ) over all directions v, i.e.
def

SI(f, κ) =supSI(f, v, κ).

v:kvk2 =1

For any constant c, we define the class of densities CSI (c, d) to consist of all d-dimensional densities f withthe property that SI(f, κ) ≤ c for all κ > 0.
Our notion of shift-invariance provides a quantitative way of capturing the intuition that the density fchanges gradually on average in every direction. Several natural classes fit nicely into this framework; forexample, we note that d-dimensional standard normal distributions are easily shown to belong to CSI (1, d).
As another example, we will show later that any d-dimensional isotropic log-concave distribution belongs
to CSI (Od (1), d).
Many distributions arising in practice have light tails, and distributions with light tails can in general belearned more efficiently. To analyze learning shift-invariant distributions in a manner that takes advantage
of light tails when they are available, while accommodating heavier tails when necessary, we define classes
with different combinations of shift-invariant and tail behavior. Given a nonincreasing function g : R+ →[0, 1] which satisfies limt→+∞ g(t) = 0, we define the class of densities CSI (c, d, g) to consist of those2f ∈ CSI (c, d) which have the additional property that for all t > 0, it holds that
Pr [||x − µ|| > t] ≤ g(t),x←fwhere µ ∈ Rd is the mean of the distribution f.
As motivation for its study, we feel that CSI (c, d, g) is a simple and easily understood class that exhibits
an attractive tradeoff between expressiveness and tractability. As we show, it is broad enough to include
distributions of central interest such as multidimensional isotropic log-concave distributions, but it is also
limited enough to admit efficient noise-tolerant density estimation algorithms.
Our density estimation framework. We recall the standard notion of density estimation with respect tototal variation distance. Given a class C of densities over Rd , a density estimation algorithm for C is given
access to i.i.d. draws from f , where f ∈ C is the unknown target density to be learned. For any f ∈ C, given
any parameter ε > 0, after making some number of draws depending on d and ε the density estimation
algorithm must output a description of a hypothesis density h over Rd which, with high probability overthe draws from f , satisfies dTV (f, h) ≤ ε. It is of interest both to bound the sample complexity of such analgorithm (the number of draws from f that it makes) and its running time.
Our learning results will hold even in a challenging model of noise-tolerant density estimation for aclass C. In this framework, the density estimation algorithm is given access to i.i.d. draws from f 0 , which
is a mixture f 0 = (1 − ε)f + εfnoise where f ∈ C and fnoise may be any density. (We will sometimes saythat such an f 0 is an ε-corrupted version of f . This model of noise is sometimes referred to as Huber’s
contamination model [Hub67].) Now the goal of the density estimation algorithm is to output a descriptionof a hypothesis density h over Rd which, with probability at least (say) 9/10 over the draws from f 0 , satisfies
dTV (f 0 , h) ≤ O(ε). This is a challenging variant of the usual density estimation framework, especially for
multidimensional density estimation. In particular, there are simple distribution learning problems (such aslearning a single Gaussian or product distribution over {0, 1}n ) which are essentially trivial in the noisefree setting, but for which computationally efficient noise-tolerant learning algorithms have proved to be asignificant challenge [DKK+ 16, DKK+ 18, SCV18].

1.1

Results

Our main positive result is a general algorithm which efficiently learns any class CSI (c, d, g) in the noisetolerant model described above. Given a constant c and a tail bound g, we show that any distribution in theclass CSI (c, d, g) can be noise-tolerantly learned to any error O(ε) with a sample complexity that depends
on c, g, ε and d. The running time of our algorithm is roughly quadratic in the sample complexity, and thed+2
sample complexity is Oc,d,g (1) · 1ε(see Theorem 29 in Section 5 for a precise statement of the exact
bound). These bounds on the number of examples and running time do not depend on which member ofCSI (c, d, g) is being learned.
Application: Learning multivariate log-concave densities. A multivariate density function f over Rd is
said to be log-concave if there is an upper semi-continuous concave function φ : Rd → [−∞, ∞) suchthat f (x) = eφ(x) for all x. Log-concave distributions arise in a range of contexts and have been wellstudied; see [CDSS13, CDSS14, ADLS17, ADK15, CDGR16, DKS16a] for work on density estimation of
univariate (discrete and continuous) log-concave distributions. In the multivariate case, [KS14] gave a sample complexity lower bound (for squared Hellinger distance) which implies that Ω(1/ε(d+1)/2 ) samples areneeded to learn d-dimensional log-concave densities to error ε. More recently, [DKS16b] established thefirst finite sample complexity upper bound for multivariate log-concave densities, by giving an algorithm
3that semi-agnostically (i.e. noise-tolerantly in a very strong sense) learns any d-dimensional log-concave
density using Õd (1/ε(d+5)/2 ) samples. The algorithm of [DKS16b] is not computationally efficient, and indeed, Diakonikolas et al. ask if there is an algorithm with running time polynomial in the sample complexity,
referring to this as “a challenging and important open question.” A subsequent (and recent) work of Carpenter et al. [CDSS18] showed that the maximum likelihood estimator (MLE) is statistically efficient (i.e.,achieves near optimal sample complexity). However, we note that the MLE is computationally inefficientand thus has no bearing on the question of finding an efficient algorithm for learning log-concave densities.
We show that multivariate log-concave densities can be learned in polynomial time as a special caseof our main algorithmic result. We establish that any d-dimensional log-concave density is Od (1)-shiftinvariant. Together with well-known tail bounds on d-dimensional log-concave densities, this easily yields
that any d-dimensional log-concave density belongs to CSI (c, d, g) where the tail bound function g is inverse exponential. Theorem 29 then immediately implies the following, answering the open question of[DKS16b]:
Theorem 1. There is an algorithm with the following property: Let f be a unknown log-concave density overRd and let f 0 be an ε-corruption of f . Given any error parameter ε > 0 and confidence parameter δ > 0
and access to independentdraws from f 0 , the algorithm with probability 1 − δ outputs a hypothesis densityR
h : Rd → R≥0 such that x∈Rd |f 0 (x) − h(x)| ≤ O(ε). The algorithm runs in time Õd (1/ε2d+2 ) · log2 (1/δ)
and uses Õd (1/εd+2 ) · log2 (1/δ) many samples.
While our sample complexity is quadratically larger than the optimal sample complexity for learning logconcave distributions (from [DKS16b]), such computational-statistical tradeoffs are in fact quite common(see, for example, the work of [BSZ15] which gives a faster algorithm for learning Gaussian mixture models
by using more samples).
A lower bound. We also prove a simple lower bound, showing that any algorithm
that learns shift-invariantdd-dimensional densities with bounded support to error ε must use Ω 1/ε examples. These densities may
be thought of as satisfying the strongest possible rate of tail decay as they have zero tail mass outside ofa bounded region (corresponding to g(t) = 0 for t larger than some absolute constant). This lower boundshows that a sample complexity of at least 1/εd is necessary even for very structured special cases of ourmultivariate density estimation problem.

1.2

Our approach

For simplicity, and because it is a key component of our general algorithm, we first describe how ouralgorithm learns an ε-error hypothesis when the target distribution belongs to CSI (c, d) and also has bounded
support: all its mass is on points in the origin-centered ball of radius 1/2.
In this special case, analyzed in Section 3, our algorithm has two conceptual stages. First, we smooth thedensity that we are to learn through convolution – this is done in a simple way by randomly perturbing eachdraw. This convolution uses a kernel that damps the contributions to the density coming from high-frequency
functions in its Fourier decomposition; intuitively, the shift-invariance of the target density ensures that the
convolved density (which is an average over small shifts of the original density) is close to the originaldensity. In the second conceptual stage, the algorithm approximates relatively few Fourier coefficients ofthe smoothed density. We show that an inverse Fourier transformation using this approximation still provides
an accurate approximation to the target density.1
1

We note that a simpler version of this approach, which only uses a smoothing kernel and does not employ Fourier analysis, can4Next, in Section 4, we consider the more general case in which the target distribution belongs to the class
CSI (c, d, g) (so at this point we are not yet in the noise-tolerant framework). Here the high-level idea of ourapproach is very straightforward: it is essentially to reduce to the simpler special case (of bounded support
and good shift-invariance in every direction) described above. (A crucial aspect of this transformation
algorithm is that it uses only a small number of draws from the original shift-invariant distribution; wereturn to this point below.) We can then use the algorithm for the special case to obtain a high-accuracy
hypothesis, and perform the inverse transformation to obtain a high-accuracy hypothesis for the original
general distribution. We remark that while the conceptual idea is thus very straightforward, there are anumber of technical challenges that must be met to implement this approach. One of these is that it isnecessary to truncate the tails of the original distribution so that an affine transformation of it will havebounded support, and doing this changes the shift-invariance of the original distribution. Another is thatthe transformation procedure only succeeds with non-negligible probability, so we must run this overall
approach multiple times and perform hypothesis selection to actually end up with a single high-accuracy
hypothesis.
In Section 5, we consider the most general case of noise-tolerant density estimation for CSI (c, d, g). Recall that in this setting the target density f 0 is some distribution which need not actually belong to CSI (c, d, g)but satisfies dTV (f 0 , f ) ≤ ε for some density f ∈ CSI (c, d, g). It turns out that this case can be handled using
essentially the same algorithm as the previous paragraph. We show that even in the noise-tolerant setting,our transformation algorithm will still successfully find a transformation as above that would succeed if thetarget density were f ∈ CSI (c, d, g) rather than f 0 . (This robustness of the transformation algorithm cruciallyrelies on the fact that it only uses a small number of draws from the given distribution to be learned.) We
then show that after transforming f 0 in this way, the original algorithm for the special case can in fact learn
the transformed version of f 0 to high accuracy; then, as in the previous paragraph, performing the inverse
transformation gives a high-accuracy hypothesis for f 0 .
In Section 6 we apply the above results to establish efficient noise-tolerant learnability of log-concave
densities over Rd . To apply our results, we need to have (i) bounds on the rate of tail decay, and (ii) shiftinvariance bounds. As noted earlier, exponential tail bounds on d-dimensional log-concave densities arewell known, so it remains to establish shift-invariance. Using basic properties of log-concave densities, inSection 6 we show that any d-dimensional isotropic log-concave density is Od (1)-shift-invariant. Armedwith this bound, by applying our noise-tolerant learning result (Theorem 29) we get that any d-dimensional
isotropic log-concave density can be noise-tolerantly learned in time Õd (1/ε2d+2 ), using Õd (1/εd+2 ) samples. Log-concave distributions are shift-invariant even if they are only approximately isotropic. We showthat general log-concave distributions may be learned by bringing them into approximately isotropic position
with a preprocessing step, borrowing techniques from [LV07].
The lower bound. As is standard, our lower bound (proved in Section 7) is obtained via Fano’s inequality.
We identify a large set F of bounded-support shift-invariant d-dimensional densities with the following two
properties: all pairs of densities from F have KL-divergence that is not too big (so that they are hard to tellapart), but also have total variation distance that is not too small (so that a successful learning algorithm isrequired to tell them apart). The members of F are obtained by choosing functions that take one of twovalues in each cell of a d-dimensional checkerboard. The two possible values are within a small constantfactor of each other, which keeps the KL divergence small. To make the total variation distance large, wechoose the values using an error-correcting code – this means that distinct members of F have different
be shown to give a similar, but quantitatively worse, results, such as a sample complexity of essentially 1/ε2d when g(t) is zero
outside of a bounded region. However, this is worse than the lower bound of Ω(1/εd ) by a quadratic factor, whereas our algorithm
essentially achieves this optimal sample complexity.

5

values on a constant fraction of the cells, which leads to large variation distance.

1.3

Related work

The most closely related work that we are aware of was mentioned above: [HK92] obtained bounds similar
to ours for using kernel methods to learn densities that belong to various Sobolev spaces. As mentioned
above, these results do not directly apply for learning densities in CSI (c, d, g) because of the possibilityof jump discontinuities. [HK92] also proved a lower bound on the sample complexity of algorithms thatcompute kernel density estimates. In contrast our lower bound holds for any density estimation algorithm,
kernel-based or otherwise.
The assumption that the target density belongs to a Besov space (see [Kle09]) makes reference to theeffect of shifts on the distribution, as does shift-invariance. We do not see any obvious containments between
classes of functions defined through shift-invariance and Besov spaces, but this is a potential topic for further
research.
Another difference with prior work is the ability of our approach to succeed in the challenging noisetolerant learning model. We are not aware of analyses for density estimation of densities belonging to
Sobolev or Besov spaces that extend to the noise-tolerant setting in which the target density is only assumed
to be close to some density in the relevant class.
As mentioned above, shift-invariance was used in the analysis of algorithms for learning discrete probability distributions in [BX99, DDO+ 13]. Likewise, both the discrete and continuous Fourier transforms
have been used in the past to learn discrete probability distributions [DKS16c, DKS16d, DDKT16].

2

Preliminaries

We write B(r) to denote the radius-r ball in Rd , i.e. B(r) = {x ∈ Rd : x21 + · · · + x2d ≤ r2 }. If f is a
probability density over Rd and S ⊂ Rd is a subset of its domain, we write fS to denote the density of f
conditioned on S.

2.1

Shift-invariance

Roughly speaking, the shift-invariance of a distribution measures how much it changes (in total variation
distance) when it is subjected to a small translation. The notion of shift-invariance has typically been used
for discrete distributions (especially in the context of proving discrete limit theorems, see e.g. [CGS11] andmany references therein). We give a natural continuous analogue of this notion below.
Definition 2. Given a probability density f over Rd , a unit vector v, and a positive real value κ, we saythat the shift-invariance of f in direction v at scale κ, denoted SI(f, v, κ), isZ
def 1
SI(f, v, κ) = · supf (x + κ0 v) − f (x) dx.
(1)κ κ0 ∈[0,κ] Rd
Intuitively, if SI(f, v, κ) = β, then for any direction (unit vector) v the variation distance between f anda shift of f by κ0 in direction v is at most κβ for all 0 ≤ κ0 ≤ κ. The factor κ1 in the definition means that
SI(f, v, κ) does not necessarily go to zero as κ gets small; the effect of shifting by κ is measured relative to
κ.
Let
def
SI(f, κ) = sup{SI(f, v, κ) : v ∈ Rd , kvk2 = 1}.
6

For any constant c we define the class of densities CSI (c, d) to consist of all d-dimensional densities f withthe property that SI(f, κ) ≤ c for all κ > 0.
We could obtain an equivalent definition if we removed the factor κ1 from the definition of SI(f, v, κ),and required that SI(f, v, κ) ≤ cκ for all κ > 0. This could of course be generalized to enforce bounds onthe modified SI(f, v, κ) that are not linear in κ. We have chosen to focus on linear bounds in this paper tohave cleaner theorems and proofs.
We include “sup” in the definition due to the fact that smaller shifts can sometimes have bigger effects.
For example, a sinusoid with period ξ is unaffected by a shift of size ξ, but profoundly affected by a shiftof size ξ/2. Because of possibilities like this, to capture the intuitive notion that “small shifts do not lead tolarge changes”, we seem to need to evaluate the worst case over shifts of at most a certain size.
As described earlier, given a nonincreasing “tail bound” function g : R+ → (0, 1) which is absolutely
continuous and satisfies limt→+∞ g(t) = 0, we further define the class of densities CSI (c, d, g) to consist of
those f ∈ CSI (c, d) which have the additional property that f has g-light tails, meaning that for all t > 0, it
holds that Prx←f [||x − µ|| > t] ≤ g(t), where µ ∈ Rd is the mean of f.
Remark 3. It will be convenient in our analysis to consider only tail bound functions g that satisfy min{r ∈R : g(r) ≤ 1/2} ≥ 1/10 (the constants 1/2 and 1/10 are arbitrary here and could be replaced by anyother absolute positive constants). This is without loss of generality, since any tail bound function g which
does not meet this criterion can simply be replaced by a weaker tail bound function g ∗ which does meet thiscriterion, and clearly if f has g-light tails then f also has g ∗ -light tails.
We will (ab)use the notation g −1 (ε) to mean inf{t : g(t) ≤ ε}.
The complexity of learning with a tail bound g will be expressed in part using
Z ∞
√
def
g( z) dz.
Ig =0

We remark that the quantity Ig is the “right” quantity in the sense that the integral Ig is finite as long as√
the density has “non-trivial
decay”. More precisely, note that by Chebyshev’s inequality, g( z) = O(z −1 ).
R
Since the integral O(z −1 )dz diverges, this means that if Ig is finite, then the density f has a decay sharper
than the trivial decay implied by Chebyshev’s inequality.

2.2

Fourier transform of high-dimensional distributions

In this subsection we gather some helpful facts from multidimensional Fourier analysis.
While it is possible to do Fourier analysis over Rd , in this paper, we will only do Fourier analysis forfunctions f ∈ L1 ([−1, 1]d ).
R
Definition 4. For any function f ∈ L1 ([−1, 1]d ), we define fb : Rd → C by fb(ξ) = x∈Rd f (x) · eπi·hξ,xi dx.
Next, we recall the following standard claims about Fourier transforms of functions, which may befound, for example, in [SS95].
R
Claim 5. For f, g ∈ L1 ([−1, 1]d ) let h(x) = y∈Rd f (y) · g(x − y)dy denote the convolution h = f ∗ g off and g. Then for any ξ ∈ Rn , we have bh(ξ) = fb(ξ) · gb(ξ).
Next, we recall Parseval’s identity on the cube.
R
Claim 6 (Parseval’s identity). For f : [−1, 1]d → R such that f ∈ L2 ([−1, 1]d ), it holds that [−1,1]d f (x)2 dx =1 P
· ξ∈Zd |fb(ξ)|2 .
2d7

The next claim says that the Fourier inversion formula can be applied to any sequence in `2 (Zd ) to obtaina function whose Fourier series is identical to the given sequence.
P
Claim 7 (Fourier inversion formula). For any g : Zd → C such that ξ∈Zd |g(ξ)2 | < ∞, the functionP1πi·hξ,xi , is well defined and satisfies b
h(x) =
h(ξ) = g(ξ) for all ξ ∈ Zd .
d d · g(ξ) · eξ∈Z 2

We will also use Young’s inequality:
Claim 8 (Young’s inequality). Let f ∈ Lp ([−1, 1]d ), g ∈ Lq ([−1, 1]d ), 1 ≤ p, q, r ≤ ∞, such that1 + 1/r = 1/p + 1/q. Then kf ∗ gkr ≤ kf kp · kgkq .

2.3

A useful mollifier

Our algorithm and its analysis require the existence of a compactly supported distribution with fast decaying
Fourier transform. Since the precise rate of decay is not very important, we use the C ∞ function b :[−1, 1] → R+ as follows:

2

− x 2
1−x
c
·
e
if |x| < 10b(x) =(2)
0
if |x| = 1.
Here c0 ≈ 1.067 is chosen so that b is a pdf; by symmetry, its mean is 0. (This function has previously beenused as a mollifier [KNW10, DKN10].) The following fact can be found in [Joh15] (while it is proved onlyfor ξ ∈ Z, it is easy to see that the same proof holds if ξ ∈ R).
√
Fact 9. For b : [−1, 1] → R+ defined in (2) and ξ ∈ Z \ {0}, we have that |bb(ξ)| ≤ e− |ξ| · |ξ|−3/4 .
QLet us now define the function bd,γ : Rd → R+ as bd,γ (x1 , . . . , xd ) = γ1d · dj=1 b(xj /γ). Combining
this definition and Fact 9, we have the following claim:√− γ·t · (γ · t)−3/4 .
Claim 10. For ξ ∈ Zd with kξk∞ ≥ t, we have |bdd,γ (ξ)| ≤ e

The next fact is immediate from (2) and the definition of bd,γ :Fact 11. kbd,γ k∞ = (c0 /γ)d and as a consequence, kbd,γ k22 ≤ (c0 /γ)2d .

3

A restricted problem: learning shift-invariant distributions with bounded
supportAs sketched in Section 1.2, we begin by presenting and analyzing a density estimation algorithm for densities
that, in addition to being shift-invariant, have support bounded in B(1/2). Our analysis also captures the fact
that, to achieve accuracy ε, an algorithm often only needs the density to be learned to have shift invarianceat a scale slightly finer than ε.
Lemma 12. There is an algorithm learn-bounded with the following property: For all constant d, forall ε, δ > 0, all 0 < κ < ε < 1/2, and all d-dimensional densities f with support in B(1/2) such thatκSI(f, κ) ≤ ε/2, given access to independent draws from f , the algorithm runs in time  2d
 
 
1
1
1 1
4d
Od 2
log
log
ε
κ
κ
κδ
8

uses

Od

1
ε2

 d
 
 
1
1
1
2d
log
log
κ
κ
κδ

R
samples, and with probability 1−δ, outputs a hypothesis h : [−1, 1]d → R+ such that x∈Rd |f (x)−h(x)| ≤
ε.
 2d

Further, given any point z ∈ [−1, 1]d , h(z) can be computed in time Od log κ(1/κ)
and satisfiesd 2dh(z) ≤ Od log κ(1/κ).
dProof. Let 0 < γ := √κd , and let us define q = f ∗ bd,γ . (Here ∗ denotes convolution and bd,γ is themollifier defined in Section 2.3.) We make a few simple observations about q:(i) Since γ ≤ 1/2, we have that q is a density supported on B(1).
(ii) Since d is a constant, a draw from bd,γ can be generated in constant time. Thus given a draw from f ,one can generate a draw from q in constant time, simply by generating a draw from bd,γ and adding itto the draw from f .
(iii) By Young’s inequality (Claim 8), we have that kqk2 ≤ kf k1 · kbd,γ k2 . Noting that f is a density andthus kf k1 = 1 and applying Fact 11, we obtain that kqk2 is finite. As a consequence, the Fourier
coefficients of q are well-defined.
Preliminary analysis. We first observe that because bd,γ is supported on [−γ, γ]d , the distribution
q may be√viewed as an average of different shifts of f where each shift is by a distance√at most γ d ≤ κ. Fix any
direction v and consider a shift of f in direction v by some distance at most γ d ≤ κ. Since κSI(f, κ) ≤ε/2, we have that the variation distance between f and this shift in direction v is at most ε/2. Averaging
over all such shifts, it follows thatdTV (q, f ) ≤ ε/2.
(3)
Next, we observe that by Claim 5, for any ξ ∈ Zd , we have qb(ξ) = fb(ξ) · bdd,γ (ξ). Since f is a pdf,+bd|f (ξ)| ≤ 1, and thus we have |bq (ξ)| ≤ |bd,γ (ξ)|. Also, for any parameter k ∈ Z , define Ck = {ξ ∈ Zd :kξk∞ = k}. Let us fix another parameter T (to be determined later). Applying Claim 10, we obtain
X
X
X X
2
2
|b
q (ξ)|2 ≤|bd|bdd,γ (ξ)| ≤d,γ (ξ)|ξ:kξk∞ >T≤X√

k>T ξ:kξk∞ =k

ξ:kξk∞ >T|Ck | · e−2· γ·k · (γ · k)−3/2 ≤

k>T

X

√

(2k + 1)d · e−2· γ·k · (γ · k)−3/2 .

k>T

√
P
, then ξ:kξk∞ >T |bq (ξ)|2 ≤ 2(2T + 1)d · e−2· γ·T ·
 
 
P
2
2
(γ · T )−3/2 . If we now set T to be 4dγ · ln2 γd + γ1 · ln2 8ε , then ξ:kξk∞ >T |bq (ξ)|2 ≤ ε8 .
2
An easy calculation shows that if T ≥ 4dγ · ln2

 
d
γ

The algorithm. We first observe that for any ξ ∈ Zd , the Fourier coefficient qb(ξ) can be estimated togood accuracy using relatively few draws from q (and hence from f , recalling (ii) above). More precisely,as an easy consequence of the definition of the Fourier transform, we have:
Observation 13. For any ξ ∈ Zd , the Fourier coefficient qb(ξ) can be estimated to within additive error of
magnitude at most η with confidence 1 − β using O(1/η 2 · log(1/β)) draws from q.
9Let us define the set Low of low-degree Fourier coefficients as Low = {ξ ∈ Zd : kξk∞ ≤ T }. Thus,|Low| ≤ (2T + 1)d . Thus, using S = O(η −2 · log(T /δ)) draws from f , by Observation 13, with probability
1 − δ, we can compute a set of values {bu(ξ)}ξ∈Low such that
For all ξ ∈ Low, |bu(ξ) − qb(ξ)| ≤ η.

(4)

Recalling (ii), the sequence {bu(ξ)}ξ∈Low can be computed in O(|S| · |Low|) time. Define ub(ξ) = 0 fordξ ∈ Z \ Low. Combining (4) with this, we getXXX|bu(ξ) − qb(ξ)|2 ≤|bu(ξ) − qb(ξ)|2 +|bu(ξ) − qb(ξ)|2ξ∈Zd

ξ6∈Low

ξ∈Low

≤

X

|b
u(ξ) − qb(ξ)|2 +ξ∈Low≤ |Low| · η 2 +

ε2
8

ε2
ε2
≤ (2T + 1)d · η 2 + .
88

2Thus, setting η as η 2 = (2T + 1)−d · ε8 , we get thatX|bu(ξ) − qb(ξ)|2 ≤

ξ∈Zd

ε2
.
4

(5)

P
u(ξ)|2 < ∞. Thus, we can apply the Fourier inversion
Note that by definition ub : Zd → C satisfies ξ∈Zd |bformula (Claim 7) to obtain a function u : [−1, 1]d → C such that
Z
X

1
ε2
|u(x) − q(x)|2 dx = d ·|bu(ξ) − qb(ξ)|2 ≤,(6)2
4 · 2d[−1,1]dd
ξ∈Z

where the first equality follows by Parseval’s identity (Claim 6). By the Cauchy-Schwarz inequality,sZZ√|u(x) − q(x)|dx ≤ 2d ·|u(x) − q(x)|2 dx.
[−1,1]d[−1,1]dR
Plugging in (6), we obtain [−1,1]d |u(x) − q(x)|dx ≤ 2ε . Let us finally define h (our final hypothesis),h : [−1, 1]d → R+ , as follows: h(x) = max{0, Re(u(x))}. Note that since q(x) is a non-negative realvalue for all x, we haveZZε|h(x) − q(x)|dx ≤|u(x) − q(x)|dx ≤ .
(7)2[−1,1]d[−1,1]dR
Finally, recalling that by (3) we have dTV (f, q) ≤ 2ε , it follows that [−1,1]d |h(x) − f (x)|dx ≤ ε.
Complexity analysis. We now analyze the time and sample complexity of this algorithm as well as thecomplexity of computing h. First of
in the value of γ and recalling that d is a all,
 observe that plugging
 2

2
constant, we get that T = 4dγ · ln2 γd + γ1 · ln2 8ε = O log κ(1/κ) . Combining this with the choice10

of η (set just above (5)), we get that the algorithm uses1
S = O 2 · log
η

= Od

1
ε2



|Low|
δ





1
= O 2 · log
η

 

T
= O

δ

(2T + 1)d · log

 

ε2

T
δ





 d
 
 
1
1
1
2d
log
log
κ
κ
κδ

draws from p. Next, as we have noted before, computing the sequence {bu(ξ)} takes time
 
 d
  
1
1
1
2d
log
log
Td
O(S · |Low|) = Odκκκδ  2d
 
 
1 1
1
1
4d
= Od 2
log
.
log
ε
κ
κ
κδ


1
ε2

 2d

To compute the function u (and hence h) at any point x ∈ [−1, 1]d takes time O(|Low|) = Od log κ(1/κ).
dThis is because the Fourier inversion formula (Claim 7) has at most O(|Low|) non-zero terms.
Finally, we prove the upper bound on h. If the training examples are x1 , ..., xS , then for any z ∈ [−1, 1]d ,we have!
S
X 1
X 1
1 X πihξ,xt iπi·hξ,zie· eπi·hξ,zih(z) ≤ |u(z)| =·ub(ξ) · e
·
=
S
2d
2d
t=1
ξ∈Low
ξ∈Low
!
2d
log (1/κ)
|Low|
,
= Od
≤
2d
κd
completing the proof.
With an eye towards our ultimate goal of obtaining noise-tolerant density estimation algorithms, the nextcorollary says that the algorithm in Lemma 12 is robust to noise. All the parameters have the same meaning
and relations as in Lemma 12.
Corollary 14. Let f 0 be a density supported in B(1/2)2 such that there is a d-dimensional density fsatisfying the following two properties: (i) f satisfies all the conditions in the hypothesis of Lemma 12, and(ii) f 0 is an ε-corrupted version of f , i.e. f 0 = (1 − ε)f + εfnoise for some density fnoise . Then given accesstofrom f 0 , the algorithm learn-bounded returns a hypothesis h : [−1, 1]d → R+ which satisfies
R samples
0x∈Rd |f (x) − h(x)| ≤ 2ε. All the other guarantees including the sample complexity and time complexity
remain the same as Lemma 12.
Proof. The proof of Lemma 12 can be broken down into two parts:• f can be approximated by q, and2

Looking ahead, while in general an “ε-noisy” version of f need not be supported in B(1/2), the reduction we employ will infact ensure that we only need to deal with noisy distributions that are in fact supported in B(1/2).

11• q can be learned.
The argument that q can be learned only used two facts about it:• it is supported in [−1, 1]d , and• it has few nonzero Fourier coefficients.
So, now consider the distribution q 0 = f 0 ∗ bd,γ where bd,γ is the same distribution as in Lemma 12. Because
q 0 is the result of convolving f 0 (a density supported in B(1/2) with bd,γ , it is supported in [−1, 1]d , andhas the same Fourier concentration property that we used for q. Thus, the algorithm will return a hypothesis
distribution h(x) such that the analogue of (7) holds, i.e.
Zε|h(x) − q 0 (x)|dx ≤ .
(8)2[−1,1]dRecalling that the density f 0 can be expressed as (1 − ε)f + εfnoise where fnoise is some density supportedin B(1/2), we now havedTV (q 0 , f 0 ) = dTV (f 0 ∗ bd,γ , f 0 ) = dTV ((1 − ε)f ∗ bd,γ + εfnoise ∗ bd,γ , (1 − ε)f + εfnoise )≤ (1 − ε)dTV (f ∗ bd,γ , f ) + εdTV (fnoise ∗ bd,γ , fnoise )≤ ε/2 + ε ≤ 3ε/2.
The penultimate inequality uses (3) and the fact that the total variation distance between any two distributions
is bounded by 1. Combining the above with (8), the corollary is proved.

4

Density estimation for densities in CSI (c, d, g)

Fix any nonincreasing tail bound function g : R+ → [0, 1] which satisfies limt→+∞ g(t) = 0 and thecondition min{r ∈ R : g(r) ≤ 1/2} ≥ 1/10} of Remark 3 and any constant c ≥ 1. In this section we prove
the following theorem which gives a density estimation algorithm for the class of distributions CSI (c, d, g):
Theorem 15. For any c, g as above and any d ≥ 1, there is an algorithm with the following property: Let f
be any target density (unknown to the algorithm) which belongs to CSI (c, d, g). Given any error parameter0 < ε < 1/2 and confidence parameter δ > 0 and access to independentRdraws from f , the algorithm withprobability 1 − O(δ) outputs a hypothesis h : [−1, 1]d → R≥0 such that x∈Rd |f (x) − h(x)| ≤ O(ε).
The algorithm runs in time!
!
 2d+2
 −1 
 −1 
1
g
(ε)
g
(ε)
1
Oc,d
(g −1 (ε))2d
log4d
log
+ Ig logεεεδδand uses
Oc,d

(g

−1

!
!
 d+2
 −1 
 −1 
1g (ε)
12d g (ε)
(ε))
log
log
+ Ig log
ε
ε
εδ
δ
d

samples.

124.1

Outline of the proof

Theorem 15 is proved by a reduction to Lemma 12. The main ingredient in the proof of Theorem 15 isa “transformation algorithm” with the following property: given as input access to i.i.d. draws from any
density f ∈ CSI (c, d, g), the algorithm constructs parameters which enable draws from the density f to betransformed into draws from another density, which we denote r. The density r is obtained by approximating
f after conditioning on a non-tail sample, and scaling the result so that it lies in a ball of radius 1/2.
Given such a transformation algorithm, the approach to learn f is clear: we first run the transformation
algorithm to get access to draws from the transformed distribution r. We then use draws from r to run thealgorithm of Lemma 12 to learn r to high accuracy. (Intuitively, the error relative to f of the final hypothesis
density is O(ε) because at most O(ε) comes from the conditioning and at most O(ε) from the algorithmof Lemma 12.) We note that while this high-level approach is conceptually straightforward, a number oftechnical complications arise; for example, our transformation algorithm only succeeds with some nonnegligible probability, so we must run the above-described combined procedure multiple times and perform
hypothesis testing to identify a successful final hypothesis from the resulting pool of candidates.
The rest of this section is organized as follows: In Section 4.2 we give various necessary technicalingredients for our transformation algorithm. We state and prove the key results about the transformation
algorithm in Section 4.3, and we use the transformation algorithm to prove Theorem 15 in Section 4.4.

4.2

Technical ingredients for the transformation algorithm

As sketched earlier, our approach will work with a density obtained by conditioning f ∈ SI(c, d) on lying
in a certain ball that has mass close to 1 under f . While we know that the original density f ∈ SI(c, d) has
good shift-invariance, we will further need the conditioned distribution to also have good shift-invariance inorder for the learn-bounded algorithm of Section 3 to work. Thus we require the following simple lemma,
which shows that conditioning a density f ∈ SI(c, d) on a region of large probability cannot hurt its shift
invariance too much.
Lemma 16. Let f ∈ SI(c, d) and let B be a ball such that Prx∼f [x ∈ B] ≥ 1 − δ where δ < 1/2. If fB
is the density of f conditioned on B, then, for all κ > 0, SI(fB , κ) ≤ 4δ
κ + 2c.
Proof. Let v be any unit vector in Rd . Note that f can be expressed as (1 − δ)fB + δ · ferr where ferr issome other density. As a consequence, for any κ > 0, using the triangle inequality we have thatZZ
|f (x) − f (x + κv)|dx ≥ (1 − δ) |fB (x) − fB (x + κv)|dxxxZ− δ |ferr (x) − ferr (x + κv)|dx.
x

Since f ∈ CSI (c, d) the left hand side is at most cκ, whereas the subtrahend on the right hand side is triviallyat most 2δ. Thus, we get
Z
2δ
cκ
+
,
(9)
|fB (x) − fB (x + κv)|dx ≤1−δ 1−δxcompleting the proof.
If f is an unknown target density then of course its mean is also unknown, and thus we will need toapproximate it using draws from f . To do this, it will be helpful to convert our condition on the tails of f tobound the variance of ||x − µ||, where x ∼ f.
13

Lemma 17. For any f ∈ CSI (c, d, g), we have Ex∼f [||x − µ||2 ] ≤ Ig .
R∞
R∞ √Proof. We have Ex∼f [||x − µ||2 ] = 0 Prx∼f [||x − µ||2 ≥ z] dz ≤ 0 g( z) dz = Ig .
The following easy proposition gives a guarantee on the quality of the empirical mean:b is its empirical estimate based on
Lemma 18. For any f ∈ CSI (c, d, g), if µ ∈ Rd is the mean of f and µ
M samples, then for any t > 0 we haveIgb ||2 ≥ t ≤Pr ||µ − µ.
Mt
Proof. If x1 , . . . , xM are independent draws from f , then
x1 + . . . + xM 2
2b || ] = E µ −E[||µ − µ
M


M
X 1
Ig
2
=
E µ − xi
=
,
2
M
M
i=1

where the last inequality is by Lemma 17. Applying Markov’s inequality on the left hand side, we get thestated claim.

4.3

Transformation algorithm

Lemma 19. There is an algorithm compute-transformation such that given access to samples from f ∈CSI (c, d, g) and an error parameter 0 < ε < 1/2, the algorithm takes O(Ig ) samples from f and withprobability at least 9/10 produces a vector µ̃ ∈ Rd and a real number t with the following properties:√
1. For Bt = {x : ||x − µ̃|| ≤ t}, we have Prx∼f [x ∈ Bt ] ≥ 1 − ε.
2. t = O(g −1 (ε)2 ),3. For all κ > 0, the density fBt satisfies SI(fBt , κ) ≤ 4ε
κ + 2c.
Proof. For M = 100Ig , the algorithm compute-transformation simply works as follows: set µ̃ to be
the empirical mean of the M samples, and t = 2((g −1 (ε))2 + 1/10). (Note that by Remark 3 we havet = Θ(g −1 (ε)2 ).). Let µ denote the true mean of f . First, by Lemma 18, with probability at least 0.9, theb will be close to the true mean µ in the following sense:
empirical mean µb ||2 ≤||µ − µ1.
10

Let us assume for the rest of the proof that this happens; fix any such outcome and denote it µ̃.
We have||x − µ̃||2 ≤ 2(||x − µ||2 + ||µ − µ̃||2 ) ≤ 2(||x − µ||2 + 1/10)

and soPrx∈f [||x − µ̃||2 > t] ≤ Prx∈f [2(||x − µ||2 + 1/10) > t] = Pr[kx − µk2 ≥ g −1 (ε)] ≤ ε.
Applying Lemma 16 completes the proof.
14

(10)

The following proposition elaborates on the properties of the output of the transformation algorithm.
Lemma 20. Let f ∈ CSI (c, d, g), ε > 0, µ̃ ∈ Rd , and t ∈ R satisfy the properties stated in Lemma 19.
Consider the density fscond defined by
√

def √
fscaled (x) = 2 t · f 2 t · (x + µ̃)
deffscond (x) = fscaled,B(1/2) (x),where fscaled,B(1/2) is the result of conditioning fscaled on membership in B(1/2). Then the density fscond (x)satisfies the following properties:
1. The density fscond is supported in the ball B(1/2).
2. For all ε < 1/2 and κ > 0, the density fscond satisfies
SI(fscond , κ) ≤√4ε
+ 4c t.
κ

Proof. First, it is easy to verify that function fscond defined above is indeed a density. Item 1 is enforced by
fiat. Now, for any direction v, we haveZ
1
fscaled (x + κ0 v) − fscaled (x) dx
SI(fscaled , v, κ) = · supκ κ0 ∈[0,κ] Rd
√
Z
√
√
2 t
=
f (2 t(x + κ0 v)) − f (2 tx) dx.
· supκ κ0 ∈[0,κ] Rd√
Using a change of variables, u = 2 tx, we get
1
SI(fscaled , v, κ) = · supκ κ0 ∈[0,κ]Z√
f (u + κ0 2 tv) − f (u) duRdZ1= · sup√f (u + κ0 v) − f (u) duκ κ0 ∈[0,2 tκ] Rd√√√= 2 t · SI(f, v, 2 tκ) ≤ 2c t.

(11)√
The last inequality uses that f ∈ CSI (c, d, g). Inequality (11) implies that fscaled ∈ CSI (2c t, d, g). Now,
Prx∼fscaled (x ∈ B(1/2)) = Prx∼f (x ∈ Bt ) ≥ 1 − ε, so applying Lemma 16 completes the proof.

4.4

Proof of Theorem 15

We are now ready to prove Theorem 15. Consider the following algorithm, which we call constructcandidates:
1. Run the transformation algorithm compute-transformation D := O(ln(1/δ)) many times (with parameter ε each time). Let (µ̃(i) , t) be the output that it produces on the i-th run, where t = O(g −1 (ε)2 ).
√
(i)
(i)
2. For each i ∈ [D], let Bt = {x : ||x − µ̃|| ≤ t} and fscond be the density defined from (µ̃(i) , t) asin Lemma 20.
15

Before describing the third step of the algorithm, we observe that given the pair (µ̃(i) , t) it is easy to(i)check whether any given x ∈ Rd belongs to Bt . We further make the following observations:(i)• If Prx∼f [x ∈ Bt ] ≥ 1/2, then with probability at least 1/2 a draw from f can be used as a draw(i)from fB (i) . In this case, via rejection sampling, it is easy to very efficiently simulate draws from fscondtgiven access to samples from f (the average slowdown is at most a factor of 2). Note that if (µ̃(i) , t)(i)satisfies the properties of Lemma 19, then Prx∼f [x ∈ Bt ] ≥ 1 − ε and we fall into this case.
(i)• On the other hand, if Prx∼f [x ∈ Bt ] < 1/2, then it may be inefficient to simulate draws from(i)fscond . But any such i will not satisfy the properties of Lemma 19, so if rejection sampling is inefficient(i)
to simulate draws from fscond then we can ignore such an i in what follows.
With this in mind, the third and fourth steps of the algorithm are as follows:(i)

3. For each i ∈ [D],3 run the algorithm learn-bounded using m samples from fscond , where m =(i)m(ε, δ, d) is the sample complexity of learn-bounded from Lemma 12. Let hscond be the resulting
hypothesis that learn-bounded outputs.
4. Finally, for each i ∈ [D] output the hypothesis obtained by inverting the mapping of Lemma 20, i.e.


1
def 1
(i)
(i)
(i)
√ · (x − µ̃ ) .
(12)h (x) = √ · hscond
2 t
2 t
Thus the output of construct-candidate is a D-tuple of hypotheses (h(1) , . . . , h(D) ).
We now analyze the construct-candidate algorithm. Given Lemma 19 and Lemma 20, it is not difficult
to show that with high probability at least one of the hypotheses that it outputs has error O(ε) with respect
to f :
R
Lemma 21. With probability at least 1 − O(δ), at least one h(i) has x |h(i) (x) − f (x)|dx ≤ O(ε).
Proof. It is immediate from Lemma 19 and the choice of D that with probability 1 − δ at least one triple(µ̃(i) , t) satisfies the properties of Lemma 19. Fix i0 to be an i for which this holds.
Given any i ∈ [D], it is easy to carry out the check for whether rejection sampling is too inefficient in(i)simulating fscond in such a way that algorithm learn-bounded will indeed be run to completion (as opposed(i0 )to being terminated) on fscond with probability at least 1 − δ, so we henceforth suppose that indeed learn(i0 )0bounded is actually run to completion on fscond . Since (µ̃(i ) , t) satisfies the properties of Lemma 19, by(i0 )Lemma 20, taking κ = min{ε/2, ε/(4g −1 (ε)c)}) the density fscond satisfies the required conditions forLemma 12 to apply with that choice of κ. The following simple proposition implies that h(i) is likewiseO(ε)-close to fBt :Proposition 22. Let f and g be two densities in Rd and let x 7→ A(x − z) be any invertible linear
transformation over Rd . Let fA (x) = det(A) · f (A(x − z)) and gA (x) = det(A) · g(A(x − z)) be thedensities from f and g under this transformation. Then dTV (f, g) = dTV (fA , gA ).
3

Actually, as described above, this and the fourth step are done only for those i for which rejection sampling is not too inefficient(i)
in simulating draws from fscond given draws from f ; for the other i’s, the run of learn-bounded is terminated.

16

Proof.
Z

Z

det(A)|f (A(x − z)) − g(A(x − z))|dx|fA (x) − gA (x)|dx =dTV (fA , gA ) =Zxx|f (z) − g(z)|dz = dTV (f, g),=zwhere the penultimate equality follows by a linear transformation of variables.
It remains only to observe that by property 1 of Lemma 19 the density fBt is ε-close to f , and then by
the triangle inequality we have that h(i) is O(ε)-close to f . This gives Lemma 21.
Tracing through the parameters, it is straightforward to verify that the sample and time complexities of
construct-candidates are as claimed in the statement of Theorem 15. These sample and time complexities dominate the sample and time complexities of the remaining portion of the algorithm, the hypothesis
selection procedure discussed below.
All that is left is to identify a good hypothesis from the pool of D candidates. This can be carriedout rather straightforwardly using well-known tools for hypothesis selection. Many variants of the basic
hypothesis selection procedure have appeared in the literature, see e.g. [Yat85, DK14, AJOS14, DDS12,
DDS15]). The following is implicit in the proof of Proposition 6 from [DDS15]:
Proposition 23. Let D be a distribution with support contained in a set W and let Dε = {Dj }M
j=1 bea collection of M hypothesis distributions over W with the property that there exists i ∈ [M ] such thatdTV (D, Di ) ≤ ε. There is an algorithm SelectD which is given ε and a confidence parameter δ, and isprovided with access to (i) a source of i.i.d. draws from D and from Di , for all i ∈ [M ]; and (ii) a (1 + β)“approximate evaluation oracle” evalDi (β), for each i ∈ [M ], which, on input w ∈ W , deterministicallyi (w)outputs D̃iβ (w) such that the value D1+β
≤ D̃iβ (w) ≤ (1 + β) · Di (w). Further, (1 + β)2 ≤ (1 + ε/8). The
SelectD algorithm has the following behavior: It makes m = O (1/ε2 ) · (log M + log(1/δ)) draws fromD and from each Di , i ∈ [M ], and O(m) calls to each oracle evalDi , i ∈ [M ]. It runs in time poly(m, M )(counting each call to an evalDi oracle and draw from a Di distribution as unit time), and with probability
1 − δ it outputs an index i? ∈ [M ] that satisfies dTV (D, Di? ) ≤ 6ε.
As suggested above, the remaining step is to apply Proposition 23 to the list of candidate hypothesis h(i)which satisfies the guarantee of Lemma 21. However, to bound the sample and time complexity of runningthe procedure Proposition 23, we need to bound the complexity both of sampling from {h(i) }i∈[D] as well asof constructing approximate evaluation oracles for these measures.4 In fact, we will first construct densities
out of the measures {h(i) }i∈[D] and show how to both efficiently sample from these measures as well as
construct approximate evaluation oracles for these densities.
(i)
Towards this, let us now define Hmax as follows: Hmax = maxi∈[D] maxz∈[−1,1]n hscond (z). From−1
Lemma 12 (recall that Lemma 12 was applied with κ = min{ε/2, ε/(4g (ε)c)}) we get that!
 −1 d−1 (ε)g (ε)
gHmax = Oc,d
log2d
.
(13)
ε
ε
We will carry out the rest of our calculations in terms of Hmax .
Note that while h(i) are forced to be non-negative and thus can be seen as measures, they need not integrate to 1 and thus neednot be densities.
4

17

R
(i)
Observation 24. For any i ∈ [D], x∈[−1,1]d hscond (x)dx can be estimated to additive accuracy ±ε and 2
confidence 1 − δ in time Od Hεmax· log(1/δ) .
2(i)

Proof. First note that it suffices to estimate the quantity Ex∈[−1,1]d [hscond (x)] to additive error ε/2d . How(i)ever, this can be estimated using the trivial random sampling algorithm. In particular, as hscond (x) ∈(i)
2 .
[0, Hmax ], the variance of the simple unbiased estimator for Ex∈[−1,1]d [hscond (x)] is also bounded by Hmax
This finishes the proof.
Note that, while the algorithm of Observation 24 does random sampling, this sampling is not from f , so
it adds nothing to the sample complexity of the learning algorithm.
R
Next, for i ∈ [D], let us define the quantity Zi to be Zi = x h(i) (x)dx. Since the functions h(i) and(i)hscond are obtained from each other by linear transformations (recall (12)), we get that that


Z
√
1
(i)
(i)
√ · (x − µ̃ ) dx.
2 tZi = hscond
2 t
x
(i)

We now define the functions H (i) and Hscond as(i)1· (x − µ̃(i) ))hscond ( 2√1h(i) (x)(i)
tand Hscond (x) =· √ .
H (x) =
Zi
Zi
2 t
(i)

(i)

Observe that the functions H (i) and Hscond are densities (i.e. they are non-negative and integrate to 1).
First, we will show that it suffices to run the procedure SelectD Ron the densities H (i) . To see this, note that
Lemma 21 says that there exists i ∈ [D] such that h(i) satisfies x |h(i) (x) − f (x)| = O(ε). For such an i,
Zi ∈ [1 − O(ε), 1 + O(ε)]. Thus, we have the following corollary.
R
Corollary 25. With probability at least 1 − δ, at least one H (i) satisfies x |H (i) (x) − f (x)| = O(ε).
Further, for such an i, Zi ∈ [1 − O(ε), 1 + O(ε)].
Thus, it suffices to run the procedure SelectD on the candidate distributions {H (i) }i∈[D] . The next
proposition shows that the densities {H (i) }i∈[D] are samplable.
Proposition 26. A draw from the density H (i) (x) can be sampled in time O(Hmax /Zi ).
(i)

(i)

Proof. First of all, note that it suffices to sample from Hscond since H (i) and Hscond are linear transforma(i)tions of each other. However, sampling from Hscond is easy using rejection sampling. More precisely, the(i)(i)
distribution Hscond is supported on [−1, 1]d . We sample from Hscond as follows:
1. Let C = [−1, 1]d × [0, Hmax ]. Sample a uniformly random point z 0 = (z1 , . . . , zd+1 ) from C.
(i)

2. If zd+1 ≤ hscond (z1 , . . . , zd ), then return the point z = (z1 , . . . , zd ).
3. Else go to Step 1 and repeat.
Now note that conditioned on returning a point in step 2, the point z is returned with probability proportional(i)(i)
to hscond (z). Thus, the distribution sampled by this procedure is indeed Hscond (z). To bound the probabilityof success, note that the total volume of C is 2d × Hmax . On the other hand, step 2 is successful only if z 0
falls in a region of volume Zi . This finishes the proof.
18

The next proposition says that if Zi ≥ 1/2, then there is an approximate evaluation oracle for the density
H (i) .
Proposition 27. Suppose Zi ≥ 1/2. Then thereis a (1+ O(ε))- approximate evaluation oracle for H (i)2which can be computed at any point w in time O Hεmax.
2
Proof. Note that we can evaluate h(i) at any point w exactly and thus the only issue is to estimate thenormalizing factor Zi . Note that since Zi ≥ 1/2 , estimating Zi to within an additive O(ε) gives us a (1 +O(ε)) multiplicativeto Zi and hence to H (i) (w) at any point w. However, by Observation 24,
 approximation
2
Hmax
this takes time O ε2 , concluding the proof.
We now apply Proposition 23 as follows.
bi .
1. For all i ∈ [D], estimate Zi using Observation 24 up to an additive error ε. Let the estimates be Z
b i ≥ 1/2}.
2. Let us define Sfeas = {i ∈ [D] : L
3. We run the routine SelectD on the densities {H (i) }i∈Sfeas . To sample from a density H (i) , we use
Proposition 26. We also construct a β = ε/32 approximation oracle for each of the densities H (i)using Proposition 27. Return the output of SelectD .
The correctness of the procedure follows quite easily.
R Namely, note that Corollary 25 implies that there is
one i such that both Zi ∈ [1 − O(ε), 1 + O(ε)] and x |H (i) (x) − f (x)| = O(ε). Thus such an i will be in
Sfeas . Thus, by the guarantee of SelectD , the output hypothesis is O(ε) close to f .
We now bound the sample complexity and time complexity of this hypothesis selection portion of thealgorithm. First of all, the number of samples required from f for running SelectD is O((1/ε2 )·(log(1/δ)+
d2 log d + log log(1/δ)) = O((1/ε2 ) · (log(1/δ) + d2 log d). This is clearly dominated by the sample
complexity of the previous parts. To bound the time complexity, note that the time complexity of invoking
the sampling oracle for any H (i) (i ∈ Sfeas ) is dominated by the time complexity of the approximate oracle
2 /ε2 . The total number of calls to the sampling as well as evaluation oracle is upper
which is 2O(d) · Hmax
bounded by ε12 (D log D + D log(1/δ)). Plugging in the value of Hmax as well as D, we see that the total
time complexity is dominated by the bound in the statement of Theorem 15. This finishes the proof.

5

Noise-tolerant density estimation for CSI (c, d, g)

Fix any nonincreasing tail bound function g : R+ → [0, 1] which satisfies limt→+∞ g(t) = 0 and thecondition of Remark 3 and any constant c ≥ 1. In this section we prove our main result, Theorem 29, which
gives a noise tolerant density estimation algorithm for CSI (c, d, g). We first recall the precise model of noise
we consider in this paper.
Definition 28. For two densities f and f 0 ∈ Rd , we say that f 0 is an ε-corruption of f if f 0 can be expressed
as f 0 = (1 − ε) · f + ε · ferr where ferr is a density in Rd .
This model of noise is sometimes referred to as Huber’s contamination model [Hub67].
Theorem 29 (Noise-tolerant density estimation for CSI (c, d, g).). For any c, g as above and any d ≥ 1,
there is an algorithm with the following property: Let f be any density (unknown to the algorithm) which
belongs to CSI (c, d, g) and let f 0 be an ε-corruption of f . Given ε and any confidence parameter δ > 0
19

and access to independent draws
from f 0 , the algorithm with probability 1 − O(δ) outputs a hypothesisRh : [−1, 1]d → R≥0 such that x∈Rd |f 0 (x) − h(x)| ≤ O(ε).
The algorithm runs in time!
!
 2d+2
 −1 
 −1 
1
g (ε)
−1
2d 1
4d g (ε)
Oc,d exp(O(Ig )) · (g (ε))
log
+ Ig log
log
ε
ε
εδ
δ
and usesOc,dexp(O(Ig )) ·

(g

−1

!
!
 d+2
 −1 
 −1 
11g (ε)
2d g (ε)(ε))
log+ Ig log
log
ε
ε
εδ
δ
d

samples.
Theorem 29 is identical to Theorem 15 except that now the target density from which draws are received
is f 0 , which is an ε-corruption of f ∈ CSI (c, d, g), rather than f itself, and the dependence on Ig is exponential. (On the other hand, if f is isotropic, then recalling Lemma 17 the function g can be taken to besuch that Ig = Ex∼f [||x − µ||2 ] = d, so that exp(O(Ig )) = O(1) and the complexity is the same as in thenoise-free case.) This requires essentially no changes in the algorithm and fortunately most of the analysis
from earlier can also be reused in a fairly black-box way. We briefly explain how the analysis of Section 4
can be augmented to handle having access to draws from f 0 rather than f .
Proof of Theorem 29. Without loss of generality, we can assume that ε ≤ 1/10, as otherwise any hypothesis
is trivially O(ε)-close to the target density. Recall that f 0 can be expressed as f 0 = (1 − ε) · f + ε · ferr
for some density ferr ∈ Rd . Since ε ≤ 1/10, this means that with probability 9/10, a random sample
from f 0 is in fact distributed exactly as a random sample from f . We now revisit the steps in the algorithm
construct-candidates (Proof of Theorem 15) and briefly sketch why sample access to f 0 instead of f
suffices.
Step 1: Note that each invocation of the algorithm compute-transformation is supposed to draw O(Ig )samples from f . We have access to f 0 rather than f , but since each sample from f 0 , with probability at least
9/10, is a sample from f , with probability at least exp(−O(Ig )), a run of compute-transformation withsamples from f 0 is the same as a run with samples from f . So now in Step 1, the algorithm computetransformation is run exp(−O(Ig )) · ln(1/δ)) many times rather than O(ln(1/δ)) many times as in theoriginal version (this accounts for the additional exp(O(Ig )) factor in the bounds of Theorem 29 versus
Theorem 15).
(i)
Step 2: This goes exactly as before with no changes. In particular, for every i, if ftrans is the true density(i)(i)
obtained by the ith transformation, then we have sample access to f 0 trans which is an ε-corruption of ftrans .
(i)
Step 3: In Step 3, we run the routine learn-bounded as usual. In particular, let us assume that ftrans satisfies the conditions of Lemma 12. Then Corollary 14, which established noise-tolerance of learn-bounded,
(i)
(i)
(i)
implies that with sample access to f 0 trans , the resulting hypothesis h0 trans is 2ε-close to ftrans .
It is easy to verify that Step 4 and the subsequent steps in hypothesis testing can go on exactly as before
with sample access to f 0 instead of f . In particular, the hypothesis testing routine SelectD will output ahypothesis h(i) which is 2ε close to f .

20

6

Efficiently learning multivariate log-concave densities

In this section we present our main application, which is an efficient algorithm for noise-tolerantly learning
d-dimensional log-concave densities. We prove the following:
Theorem 30 (Restatement of Theorem 1). There is an algorithm with the following property: Let f be aunknown log-concave density over Rd and let f 0 be an ε-corruption of f . Given any error parameter ε > 0
and confidence parameter δ > 0 and access to independent draws
from f 0 , the algorithm with probabilityR
1 − δ outputs a hypothesis density h : Rd → R≥0 such that x∈Rd |f 0 (x) − h(x)| ≤ O(ε). The algorithm
runs in time!
 2d+2
 
 
1
1
1
1
Od
log
log
log7d
ε
ε
εδ
δ
and usesOd!
 d+2
 
 
1
1
1
4d 1
log
log
log
ε
ε
εδ
δ

samples.
We will establish Theorem 30 in two stages. First, we will show that any log-concave f that is nearly
isotropic in fact belongs to a suitable class CSI (c, d); given this, the theorem follows immediately from
Theorem 29 and a straightforward tracing through of the resulting time and sample complexity bounds.
Then, we will reduce to the near-isotropic case, similarly to what was done in [LV07, BL13].
First, let us state the theorem for the well-conditioned case. For this, the following definitions will be
helpful.
Definition 31. Let Σ and Σ̃ be two positive semidefinite matrices. We say that Σ and Σ̃ are C-approximations
e 6= 0, we haveof each other (denoted by Σ ≈C Σ̃) if for every x ∈ Rn such that xT Σx
1
xT Σx
≤ C.
≤
e
C
xT Σx
Definition 32. Say that the probability distribution is C-nearly-isotropic if its covariance matrix C-approximatesI, the d-by-d identity matrix.
Theorem 33. There is an algorithm with the following property: Let f be a unknown C-nearly-isotropic
log-concave density over Rd and let f 0 be an ε-corruption of f , where C and d are constants.
Given any error parameter ε > 0 and confidence parameter δ > 0 and access to independent
R draws fromf 0 , the algorithm with probability 1 − δ outputs a hypothesis density h : Rd → R≥0 such that x∈Rd |f 0 (x) −h(x)| ≤ O(ε). The algorithm runs in time!
 2d+2
 
 
1
1
1
1
OC,d
log7d
log
log
ε
ε
εδ
δ
and usesOC,d!
 d+2
 
 
1
1
1
4d 1
log
log
log
ε
ε
εδ
δ

samples.
21

By Theorem 29, Theorem 33 is an immediate consequence of the following theorem on the shiftinvariance of near-isotropic log-concave distributions.
Theorem 34. Let f be a C-nearly-isotropic log-concave density in Rd , for constants C and d. Then, forg(t) = e−Ω(t) , there is a constant c1 = OC,d (1) such that f ∈ CSI (c1 , d, g).
Proof. The fact that f has e−Ω(t) -light tails directly follows from Lemma 5.17 of [LV07], so it remains toprove that there is a constant c1 such that f ∈ CSI (c1 , d). Because membership in CSI (c1 , d) requires that a
condition be satisfied for all directions v, rotating a distribution does not affect its membership in CSI (c1 , d).
Choose a unit vector v and κ > 0. By rotating the distribution if necessary, we may assume that v = e1 ,and our goal of showing that SI(f, e1 , κ) ≤ c1 is equivalent to showing thatZ|f (x) − f (x + κ0 e1 )|dx ≤ c1 κ
(14)for all κ0 ≤ κ.
def
We bound the integral of the LHS as follows. Fix some value of x0 = (x2 , . . . , xd ). Let us define
def

Lx0 = {(x1 , x2 , . . . , xd ) : x1 ∈ R} to be the line through (0, x2 , . . . , xd ) and (1, x2 , . . . , xd ). Since therestriction of a concave function to a line is concave, the restriction of a log-concave distribution to a line is
log-concave. SinceZ
Z Z
0|f (x) − f (x + κ e1 )| dx =|f (x1 , x2 , ..., xd ) − f (x1 + κ0 , x2 , ..., xd )| dx1 dx0
(15)
x0

x1

we are led to examine the one-dimensional log-concave measure f (·, x2 , ..., xd ). The following will be
useful for that.
Claim 35. Let ` : R → R be a log-concave measure. Then,
Z|`(t) − `(t + h)|dt ≤ 3h · max `(t).
t∈R

Proof. Log-concave measures are unimodal (see [Ibr56]). Let z be the mode of `, so that ` is non-decreasing
on the interval [−∞, z] and non-increasing in [z, ∞]. We haveZ|`(t) − `(t + h)| dt
Z z−h
Z z
Z ∞
=
|`(t) − `(t + h)| dt +|`(t) − `(t + h)| dt +|`(t) − `(t + h)| dt
−∞

z−h

Z z−h
`(t + h) − `(t) dt +

=

z

Z ∞

Z z

−∞

|`(t) − `(t + h)| dt +z−h`(t) − `(t + h) dt
z(since z is the mode of `)
Z z

Z z

=

|`(t) − `(t + h)| dt +`(t) dt +
z−h

Z z+h

z−h

`(t) dtz≤ 3h max `(t).
t∈R

22

Returning to the proof of Theorem 34, applying Claim 35 with (15), we get

Z
Z 
0
0
0
max f (x1 , x ) dx0 .
|f (x) − f (x + κ e1 )| dx ≤ 3κx0

x1 ∈Lx0(16)Now, since an isotropic log-concave distribution g satisfies g(x) ≤ K exp(−kxk) for an absolute constant K (see Theorem 5.1 of [SW14]) , our C-nearly-isotropic log-concave distribution f satisfies f (x) ≤C d K exp(−kxk) = OC,d (exp(−kxk)). Plugging this into (16), we get

Z
Z 
0
0
0
max exp(−k(x1 , x )k) dx0|f (x) − f (x + κ e1 )| dx ≤ OC,d (κ )x0 x1 ∈Lx0
Z
exp(−kx0 k) dx0 .
≤ OC,d (κ0 )x0

Since the integral converges, this finishes the proof.
Finally, we turn to the problem of learning log-concave distributions that are not C-nearly-isotropic.
To learn in this case, we need to rescale the axes if necessary to transform the distribution so that it is Cnearly-isotropic. To compute this rescaling, we will first assume that there is no noise present in the samples.
(The trick to handle noise is simple and will be discussed later.) First, as has often been observed, we may
assume without loss of generality that the covariance matrix of the target has full rank, since, otherwise thealgorithm can efficiently find the affine span of the entire distribution (possibly up to a negligible amount ofprobability mass), and the algorithm can be carried out within that lower dimensional subspace. To bring
the distribution to nearly isotropic position, we will be using ideas from [LV07]. (We require the additional
analysis below, rather than invoking their results as a black box, to cope with the fact that the mean is
unknown.)
Our starting point is the following lemma due to Lovász and Vempala [LV07].
Lemma 36. Let f be a zero-mean log-concave density on Rd . For m = O(d log3 d), if Σ denotes theb is the empirical covariance matrix from m samples of f , then,(population) covariance matrix of f and Σbwith probability 9/10, Σ is a 11/10 approximation to Σ.
Lemma 36 enables us to estimate the covariance matrix if we know the mean. To apply it when we donot, we appear to need an estimate of the mean that is especially good in directions with low variance. The
following is aimed at obtaining such an estimate.
Recall that, for a set A of real-valued functions on a common domain X, the pseudo-dimension of A,
which is denoted by Pdim(A), is the VC-dimension of the set of indicator functions, one for each a ∈ A,
for whether (x, y) satisfies a(x) ≥ y. We will use the following standard VC bound.
Lemma 37 ([Tal94]). For any set A of functions with a common domain X and ranges contained in
2[−M, M ], for any distribution D, m = O( M (Pdim(A)+log(1/γ))) suffices for a set S of m examples drawn
ε2according to D, with probability 1 − γ, to have all a ∈ A have
Ex∼D (a(x)) −1 Xa(x0 ) ≤ ε.
m 0
x ∈S

The proof of the following lemma follows a similar lemma in [KLS09].

23

Lemma 38. Fix a function b from Rd to R+ . Define au = b(u) · (u · x). The pseudo-dimension of{au : u ∈ B(1)} is O(d).
Proof. Any (x, y) satisfies au (x) ≥ y iffb(u)(u · x) ≥ y.
Thus, the set of indicator functions for au (x) ≥ y can be embedded into the set of homogeneous halfspaces
over Rd+1 , which is known to have VC-dimension O(d).
Now we are ready for the result we require on estimating the mean:2
Lemma 39. Fix any log-concave distribution f over Rd and any α > 0. For m = O d logα(d/α), with2
probability at least 3/4, a multiset S of m samples drawn i.i.d. from f satisfies, for all unit length u,| Ex∼S (u · x) − Ex∼f (u · x)|p≤ α.
Varx∼f (u · x)

(17)

Proof. Translating the distribution f translates both Ex∼S (u · x) and Ex∼f (u · x) the same way, and does
not affect Varx∼f (u · x), so we may assume without loss of generality that f has zero mean.
p Let fB be the
distribution obtained from f by conditioning the choice of x on the event that |u · x| ≤ Varx∼f (u · x) ·ln(8m)
for all unit length u, where c is a large constant. Lemma 5.17 of [LV07] implies that, for large enoughcc, the total variation distance between f and fB is 1/(8m), so that the total variation distance between mdraws from f and m draws from fB is at most 1/8. We henceforth assume that the m draws from f are infact drawn from fB , and proceed to analyze fB .
For any unit length u, define au by au (x) = √ |u·x|. Lemma 38 implies that {au : u ∈ B(1)}Varx∼f (u·x)

has pseudo-dimension O(d). Furthermore, when x is chosen from the support of fB , each au takes values
in an interval of size O(log m). Thus we may apply Lemma 37 to obtain Lemma 39.
Now we are ready to present and analyze the transformation.
Lemma 40. There is an algorithm rescale such that given access to samples from a log-concave distribution f , and an error parameter ε > 0, the algorithm takes O(d log3 d) samples from f and with probability
at least 1/2 produces a non-singular positive definite matrix Σ̃ ∈ Rd×d such that, if Σ is the covariancematrix of f , for any unit vector v,
1v T Σv≤≤ 2.
2
v T Σ̃v
Proof. For a large constant C and M = Cd log3 d the algorithm rescale first uses M examples to construct
an estimate µ̃ of the mean of f , and then uses µ̃ to use the examples estimate the covariance matrix.
Lemma 39 implies that, if C is large enough, then with probability 3/4, for all unit length v, we have|µ · v − µ̃ · v|
1
p
≤ .
10
Varx∼f [v · x]
Lemma 36 implies that, with probability 3/4 over a random i.i.d. draw of x1 , . . . , xM ∼ f, we have
1 PM
2i=1 (v · xi − v · µ)M
9/10 ≤≤ 11/10.
Varx∼f [v · x]
24

(18)

(19)

We henceforth assume that both (18) and (19) hold (this happens with probability at least 1/2), and we letµ̃ and x1 , . . . , xM denote the corresponding outcomes.
Let Σ be the true co-variance of f , and let Σ̃ be the estimate that was used (which depends on µ̃).
We have that
1 PM
2v T Σ̃v
i=1 (v · xi − v · µ̃)M=v T Σvv T Σv
1 PM(v · xi − v · µ̃)2
= M i=1
Varx∼f [v · x]PM
1(v · xi − v · µ + v · µ − v · µ̃)2
= M i=1
Varx∼f [v · x]2|v·µ−v·µ̃| PM
1 PM
1 PM
2
2
i=1 (v · xi − v · µ)i=1 |v · xi − v · µ|i=1 (v · µ − v · µ̃)
M
M
M
≤
+
+
Varx∼f [v · x]
Varx∼f [v · x]
Varx∼f [v · x]PMM
1(v · xi − v · µ)2
2|v · µ − v · µ̃| X |v · xi − v · µ|(v · µ − v · µ̃)2p
= M i=1+ p+
Varx∼f [v · x]
Varx∼f [v · x]
M Varx∼f [v · x] i=1 Varx∼f [v · x]
v
u
M
1 PM
2
u X
(v
·
x
−
v
·
µ)
2|v
·
µ
−
v
·
µ̃|
(v · xi − v · µ)2 (v · µ − v · µ̃)2
i
i=1
M
t
+ p
× M
+
≤
Varx∼f [v · x]
Varx∼f [v · x]
Varx∼f [v · x]
M Varx∼f [v · x]
i=1
v
u
M
1 PM
2
u1 X
(v
·
x
−
v
·
µ)
2|v
·
µ
−
v
·
µ̃|
(v · xi − v · µ)2 (v · µ − v · µ̃)2
i
i=1
M
t
=
+p
×
+
Varx∼f [v · x]MVarx∼f [v · x]
Varx∼f [v · x]
Varx∼f [v · x]i=1p
≤ 11/10 + 2(1/10) 11/10 + 1/100 ≤ 2,
where the second inequality is by Cauchy-Schwarz and the third inequality is by (19) and (18). Similarly,
1v T Σ̃v
M≥v T ΣvPM

i=1 (v · xi − v · µ)22|v·µ−v·µ̃| PM
i=1 |v · xi − v · µ|M−Varx∼f [v · x]p
≥ 9/10 − 2(1/10) 11/10 ≥ 1/2,

Varx∼f [v · x]completing the proof.

Proof of Theorem 30. The basic algorithm (for the noise-free setting) applies the procedure rescale from
Lemma 40 to find an estimate of the covariance matrix of f , rescales the axes so that the transformed
distribution is 2-nearly-isotropic, learns the transformed distribution, and then rescales the axes again torestore their original scales.
In the presence of noise, Lemma 40 succeeds with probability 1/2, if all the examples are not noisy. But
since the noise rate is at most ε, and we may assume without loss of generality that ε < 1/10, since thenumber of examples required in Lemma 40 is independent of ε, any invocation of the method succeeds inthe presence of noise with probability Ωd (1), which is at least some positive constant (since d is a constant).
Thus, if an algorithm performs Od (log(1/δ)) many repetitions, with probability at least 1 − δ/2 one of themwill succeed. It can therefore call the algorithm of Theorem 29 O(log(1/δ)) times, and then applying the
hypothesis testing procedure of Proposition 23 to the results, to achieve the claimed result.
25

7

Learning shift-invariant densities over Rd with bounded support requiresΩ(1/εd ) samples

In this section we give a simple lower bound which shows that Ω(1/εd ) samples are required for ε-accurate
density estimation even of shift-invariant d-dimensional densities with bounded support. As discussed in
the introduction, densities with bounded support may be viewed as satisfying the strongest possible rate of
tail decay as they have zero tail mass outside of a bounded region.
√
Theorem 41. Given d ≥ 1, there is a constant cd = Θ( d) such that the following holds: For all sufficiently
small ε, let A be an algorithm with the following property: given access to m i.i.d. samples from an arbitrary(and unknown) finitely supported density f ∈ CSI (cd , d), with probability at least 99/100, A outputs ahypothesis density h such that dTV (f, h) ≤ ε. Then m ≥ Ω((1/ε)d ).
Since an algorithm that achieves a small error with high probability can be used to achieve small error
in expectation, to prove Theorem 41 it suffices to show that any algorithm that achieves expected error O(ε)
must use Ω((1/ε)d ) samples. To establish this we use Lemma 42 (given below), which provides a lower
bound on the number of examples needed for small expected error.
To obtain the desired lower bound from Lemma 42, we establish the existence of a family F of densities
F = {f1 , . . . , fN } ∈ CSI (cd , d), where N = exp(Ω((1/ε)d )). These densities will be shown to satisfy thefollowing two properties: for any i 6= j ∈ [N ] we have (1) dTV (fi , fj ) = Ω(ε), and (2) the Kullback-Leibler
divergence DKL (fi ||fj ) is at most O(1), yielding Theorem 41.

7.1

Fano’s inequality

The main tool we use for our lower bound is Fano’s inequality, or more precisely, the following extensionof it given by [IH81] and [AB83]:
Theorem 42 (Generalization of Fano’s Inequality.). Let f1 , . . . , fN +1 be a collection of N + 1 distributions
such that for any i 6= j ∈ [N + 1], we have (i) dTV (fi , fj ) ≥ α/2, and (ii) DKL (fi ||fj ) ≤ β, where DKL
denotes Kullback-Leibler divergence. Then for any algorithm that makes m draws from an unknown target
distribution fi , i ∈ [N + 1], and outputs a hypothesis distribution f˜, there is some i ∈ [t + 1] such that if
the target distribution is fi , then


α
mβ
+
ln
2
E[dTV (f, f˜)] ≥
1−
.
2
ln N
In particular,to achieve expected error at most α/4, any learning algorithm must have sample complexity
ln N
m=Ω β .

7.2

The family of densities we analyze

Let T be a positive integer that is T = dC/εe for a large constant C. We consider probability densities over
Rd which (i) are supported on [−T, T )d , and (ii) are piecewise constant on each of the (2T )d many disjoint
unit cubes whose union is [−T, T )d . Writing A to denote the set {−T, −T + 1, . . . , −1, 0, 1, . . . , T −1}d , each of the (2T )d many disjoint unit cubes mentioned above is indexed by a unique element a =(a1 , . . . , ad ) ∈ A in the obvious way. We write cube(a) to denote the unit cube indexed by a. Givenx ∈ [−T, T )d we write a(x) to denote the unique element a ∈ A such that x ∈ cube(a).
26For any z ∈ {0, 1}A , we define a probability density fz over [−T, T )d as fz (x) = (T + za(x) )/Z, whereZ = Θ((2T )(d+1) ) is a normalizing factor so that fz is indeed a density (i.e. it integrates to 1).
It is well known (via an elementary probabilistic argument) that there is a subset S ⊂ {0, 1}A of size2Θ(|A|) such that any two distinct strings z, z 0 ∈ S differ in Θ(|A|) many coordinates. We define the set F
of densities to be F = {fz : z ∈ S}.

7.3

Membership in CSI (cd , d)

It is obvious that every density in F is finitely supported. In this subsection we prove that F ⊆ CSI (cd , d).
First, we bound the variation distance incurred by shifting along a coordinate axis:
R
Lemma 43. For any f ∈ F, i ∈ [d], and κ ∈ (0, 1), we have |f (x + κei ) − f (x)| dx ≤ Θ(κε).
Proof. We haveZZ|f (x + κei ) − f (x)| dx =Z|f (x + κei ) − f (x)| dx|f (x + κei ) − f (x)| dx +{x:xi >T −κ}{x:xi <−T }Z|f (x + κei ) − f (x)| dx.

+{x:xi ∈[−T,T −κ]}If xi < −T − κ, then f (x + κei ) = f (x) = 0. When −T − κ ≤ xi < −T , we have |f (x + κei ) − f (x)| ≤(T + 1)/Z. ThusZ|f (x + κei ) − f (x)| dx ≤ (κ(2T )d−1 )(T + 1)/Z = Θ(κε).
{x:xi <−T }A similar argument gives thatR

{x:xi >T −κ} |f (x + κei ) − f (x)| dx ≤ Θ(κε). Finally,ZZ|f (x + κei ) − f (x)| dx ={x:xi ∈[−T,T −κ]}|f (x + κei ) − f (x)| dx{x:xi ∈[−T,T −κ],dxi e−xi ≤κ}Z≤(1/Z) dx.
{x:xi ∈[−T,T −κ],dxi e−xi ≤κ}The set {x : xi ∈R[−T, T − κ], dxi e − xi ≤ κ} is made up of 2T − 1 “slabs” that are each of width κ,and consequently {x:xi ∈[−T,T −κ],dxi e−xi ≤κ} ≤ (2T − 1)κ(2T )d−1 /Z = Θ(κ/T )= Θ(κε), recalling thatZ = Θ((2T )d+1 ). This completes the proof.
Given Lemma 43, it is easy to bound the variation distance incurred by shifting in an arbitrary direction:√
R
Lemma 44. For any f ∈ F, for any unit vector v, for any κ < 1, we have |f (x+κv)−f (x)| dx ≤ cκε dfor a universal constant c > 0.

27

Proof. Writing the unit vector v asPd

j=1 vj ej , we haveZ|f (x + κv) − f (x)| dx =

Z X
d


f x + κ

i=1

≤





vj ej  − f x + κ

j=1



d Z
X

i
X

f x + κ

i=1

i
X

i−1
X


vj ej  dx

j=1





vj ej  − f x + κ

j=1

i−1
X


vj ej  dx

j=1

(triangle inequality)
=

d Z
X

|f (x + κ|vi |ei ) − f (x) | dx

(variable substitution)

i=1

≤ cκε

d
X

√
|vi | ≤ cκε d,(Lemma 43 and Cauchy-Schwarz)

i=1

completing the proof.
As an easy consequence of Lemma 44, Lemma 17 and the definition of SI(f, v, κ) we obtain the following:√
Corollary 45. There is a constant cd = Θ( d) such that for any f ∈ F and any unit vector v ∈ Rd , wehave SI(f, v) ≤ cd . Hence F ⊆ CSI (cd , d).

7.4

The upper bound on KL divergence and lower bound on variation distance

Recall that if f and g are probability density functions supported on a set S ⊆ Rd , then the Kullback LeiblerR(x)divergence between f and g is defined as DKL (f ||g) = S f (x) ln fg(x)dx. As an immediate consequenceof this definition, we have the following claim:
Claim 46. Let f, g be two densities such that for some absolute constant C > 1 we have that every x
satisfies C1 f (x) ≤ g(x) ≤ Cf (x). Then DKL (f ||g) ≤ O(1).
It is easy to see that any fi , fj in the family F of densities described in Section 7.2 are such that1f
C i (x) ≤ fj (x) ≤ Cfi (x). Thus we have:
Lemma 47. DKL (fi ||fj ) ≤ O(1) for all i 6= j ∈ [N ].
Finally, we need a lower bound on the total variation distance between any pair of elements of F:
Lemma 48. There is an absolute constant c > 0 such that, for any fu , fv ∈ F, dTV (fu , fv ) = Ω(ε).
Proof. We have dTV (fu , fv ) = (1/Z)|{a ∈ A : ua 6= va }| = (1/Z)Ω(|A|) = Ω(1/T ) = Ω(ε).

7.5

Putting it together

By Lemma 48, each pair of elements of F are separated by Ω(ε) in total variation distance. Since Lemma 47
implies that each pair of elements of F has KL-divergence O(1), Theorem 42 implies that Ω(ln |F|) =Ω((1/ε)d ) examples are needed to√achieve expected error at most O(ε). Since Corollary 45 gives that F ⊆CSI (cd , d) for a constant cd = Θ( d), this proves Theorem 41.
28

References[AB83] P. Assouad and C. Birge. Deux remarques sur l’estimation. C. R. Acad. Sci. Paris Sér. I,
296:1021–1024, 1983.
[ADK15] Jayadev Acharya, Constantinos Daskalakis, and Gautam Kamath. Optimal testing for propertiesof distributions. In Advances in Neural Information Processing Systems 28 (NIPS), pages 3591–
3599, 2015.
[ADLS17] Jayadev Acharya, Ilias Diakonikolas, Jerry Li, and Ludwig Schmidt. Sample-optimal density
estimation in nearly-linear time. In Proceedings of the Twenty-Eighth Annual ACM-SIAM Symposium on Discrete Algorithms, pages 1278–1289. SIAM, 2017.
[AJOS14] J. Acharya, A. Jafarpour, A. Orlitsky, and A.T. Suresh. Near-optimal-sample estimators forspherical gaussian mixtures. Technical Report http://arxiv.org/abs/1402.4746, ArXiV, 19 Feb
2014.
[BC91] Andrew R Barron and Thomas M Cover. Minimum complexity density estimation. IEEE
transactions on information theory, 37(4):1034–1054, 1991.
[Bes59] OV Besov. On a family of function spaces. embedding and extension theorems. Dokl. Akad.
Nauk SSSR, 126(6):1163–1165, 1959.
[BL13] Maria-Florina Balcan and Philip M Long. Active and passive learning of linear separators underlog-concave distributions. In Conference on Learning Theory, pages 288–316, 2013.
[BSZ15] Aditya Bhaskara, Ananda Suresh, and Morteza Zadimoghaddam. Sparse solutions to nonnegative linear systems and applications. In Artificial Intelligence and Statistics, pages 83–92,
2015.
[BX99] Andrew D Barbour and Aihua Xia. Poisson perturbations. ESAIM: Probability and Statistics,
3:131–150, 1999.
[CDGR16] Clément L. Canonne, Ilias Diakonikolas, Themis Gouleakis, and Ronitt Rubinfeld. Testing
shape restrictions of discrete distributions. In 33rd Symposium on Theoretical Aspects of Computer Science, STACS 2016, pages 25:1–25:14, 2016.
[CDSS13] S. Chan, I. Diakonikolas, R. Servedio, and X. Sun. Learning mixtures of structured distributions
over discrete domains. In SODA, pages 1380–1394, 2013.
[CDSS14] S. Chan, I. Diakonikolas, R. Servedio, and X. Sun. Efficient density estimation via piecewise
polynomial approximation. In STOC, pages 604–613, 2014.
[CDSS18] Timothy Carpenter, Ilias Diakonikolas, Anastasios Sidiropoulos, and Alistair Stewart. NearOptimal Sample Complexity Bounds for Maximum Likelihood Estimation of Multivariate Logconcave Densities. CoRR, abs/1802.10575, 2018.
[CGS11] L. Chen, L. Goldstein, and Q.-M. Shao. Normal Approximation by Stein’s Method. Springer,
2011.

29

[Das99] S. Dasgupta. Learning mixtures of Gaussians. In Proceedings of the 40th Annual Symposium
on Foundations of Computer Science, pages 634–644, 1999.
[DDKT16] Constantinos Daskalakis, Anindya De, Gautam Kamath, and Christos Tzamos. A size-free CLT
for poisson multinomials and its applications. In Proceedings of the 48th Annual ACM SIGACT
Symposium on Theory of Computing, STOC 2016, pages 1074–1086, 2016.
[DDO+ 13] Constantinos Daskalakis, Ilias Diakonikolas, Ryan O’Donnell, Rocco A. Servedio, and Li-Yang
Tan. Learning sums of independent integer random variables. In Foundations of Computer
Science (FOCS), 2013 IEEE 54th Annual Symposium on, pages 217–226. IEEE, 2013.
[DDS12] C. Daskalakis, I. Diakonikolas, and R.A. Servedio. Learning Poisson Binomial Distributions.
In STOC, pages 709–728, 2012.
[DDS15] A. De, I. Diakonikolas, and R. Servedio. Learning from Satisfying Assignments. In Proc.
ACM-SIAM Symposium on Discrete Algorithms (SODA), pages 478–497, 2015.
[DG85] L. Devroye and L. Györfi. Nonparametric Density Estimation: The L1 View. John Wiley &Sons, 1985.
[DK14] C. Daskalakis and G. Kamath. Faster and sample near-optimal algorithms for proper learning
mixtures of gaussians. In COLT, pages 1183–1213, 2014.
[DKK+ 16] Ilias Diakonikolas, Gautam Kamath, Daniel M. Kane, Jerry Li, Ankur Moitra, and Alistair
Stewart. Robust estimators in high dimensions without the computational intractability. In
IEEE 57th Annual Symposium on Foundations of Computer Science, FOCS, pages 655–664,
2016.
[DKK+ 18] Ilias Diakonikolas, Gautam Kamath, Daniel M. Kane, Jerry Li, Ankur Moitra, and Alistair
Stewart. Robustly learning a gaussian: Getting optimal error, efficiently. In Proceedings of
the Twenty-Ninth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA, pages 2683–
2702, 2018.
[DKN10] Ilias Diakonikolas, Daniel M Kane, and Jelani Nelson. Bounded independence fools degree-2
threshold functions. In FOCS, 2010.
[DKS16a] Ilias Diakonikolas, Daniel M. Kane, and Alistair Stewart. Efficient robust proper learning of
log-concave distributions. CoRR, abs/1606.03077, 2016.
[DKS16b] Ilias Diakonikolas, Daniel M Kane, and Alistair Stewart. Learning multivariate log-concave
distributions. In Conference on Learning Theory (COLT), pages 711–727, 2016.
[DKS16c] Ilias Diakonikolas, Daniel M Kane, and Alistair Stewart. Optimal learning via the Fouriertransform for sums of independent integer random variables. In Conference on Learning Theory, pages 831–849, 2016.
[DKS16d] Ilias Diakonikolas, Daniel M Kane, and Alistair Stewart. The Fourier transform of Poisson
multinomial distributions and its algorithmic applications. In Proceedings of the forty-eighth
annual ACM symposium on Theory of Computing, pages 1060–1073. ACM, 2016.

30

[DL12] Luc Devroye and Gábor Lugosi. Combinatorial methods in density estimation. Springer Science& Business Media, 2012.
[DLS18] A. De, P. M. Long, and R. A. Servedio. Learning sums of independent random variables withsparse collective support. FOCS, 2018.
[DS93] Ronald A DeVore and Robert C Sharpley. Besov spaces on domains in <d . Transactions of theAmerican Mathematical Society, 335(2):843–864, 1993.
[HK92] Lasse Holmström and Jussi Klemelä. Asymptotic bounds for the expected l1 error of a multivariate kernel density estimator. Journal of multivariate analysis, 42(2):245–266, 1992.
[Hub67] P. Huber. The behavior of maximum likelihood estimates under nonstandard conditions. In Proceedings of the fifth Berkeley symposium on mathematical statistics and probability, volume 1,
pages 221–233. Berkeley, CA, 1967.
[Ibr56] Ildar Abdullovich Ibragimov. On the composition of unimodal distributions. Theory of Probability & Its Applications, 1(2):255–260, 1956.
[IH81] I.A. Ibragimov and R.Z. Has’minskii. Statistical estimation, asymptotic theory (Applications ofMathematics, vol. 16). Springer-Verlag, New York, 1981.
[Joh15] S. Johnson. Saddle-point integration of c-infinity “bump” functions.
arXiv:1508.04376, 2015.

arXiv preprint[Kle09] Jussi Klemelä. Smoothing of Multivariate Data: Density Estimation and Visualization. Wiley
Publishing, 2009.
[KLS09] A. Klivans, P. Long, and R. Servedio. Learning Halfspaces with Malicious Noise. Journal of
Machine Learning Research, 10:2715–2740, 2009.
[KMV10] A. T. Kalai, A. Moitra, and G. Valiant. Efficiently learning mixtures of two Gaussians. In
STOC, pages 553–562, 2010.
[KNW10] Daniel M Kane, Jelani Nelson, and David P Woodruff. On the exact space complexity ofsketching and streaming small norms. In SODA, 2010.
[KS14] A.K.H. Kim and R.J. Samworth. Global rates of convergence in log-concave density estimation.
Available at http://arxiv.org/abs/1404.2298, 2014.
[LV07] László Lovász and Santosh Vempala. The geometry of logconcave functions and samplingalgorithms. Random Struct. Algorithms, 30(3):307–358, 2007.
[Mas97] Elias Masry. Multivariate probability density estimation by wavelet methods: Strong consistency and rates for stationary time series. Stochastic processes and their applications,
67(2):177–193, 1997.
[MV10] A. Moitra and G. Valiant. Settling the polynomial learnability of mixtures of Gaussians. In
FOCS, pages 93–102, 2010.

31

[SCV18] Jacob Steinhardt, Moses Charikar, and Gregory Valiant. Resilience: A criterion for learning in
the presence of arbitrary outliers. In 9th Innovations in Theoretical Computer Science Conference, ITCS, pages 45:1–45:21, 2018.
[Sob63] Sergej Lvovich Sobolev. On a theorem of functional analysis. Am. Math. Soc. Transl., 34:39–
68, 1963.
[SS95] Winthrop W Smith and Joanne M Smith. Handbook of real-time fast fourier transforms. IEEE,
New York, 1995.
[SW14] A. Saumard and J.A. Wellner. Log-concavity and strong log-concavity: a review. Technical
Report https://arxiv.org/pdf/1404.5886.pdf, 23 April 2014.
[Tal94] M. Talagrand. Sharper bounds for Gaussian and empirical processes. Annals of Probability,
22:28–76, 1994.
[WN07] Rebecca M Willett and Robert D Nowak. Multiscale poisson intensity and density estimation.
IEEE Transactions on Information Theory, 53(9):3171–3187, 2007.
[Yat85] Y. G. Yatracos. Rates of convergence of minimum distance estimators and Kolmogorov’s entropy. Annals of Statistics, 13:768–774, 1985.

32