An Efficient and Continuous Voronoi Density Estimator

Giovanni Luca Marchetti Vladislav Polianskii Anastasiia Varava

Florian T. Pokorny Danica Kragic

School of Electrical Engineering and Computer Science, KTH Royal Institute of
Technology Stockholm, Sweden

###### Abstract

We introduce a non-parametric density estimator deemed Radial Voronoi Density
Estimator (RVDE). RVDE is grounded in the geometry of Voronoi tessellations
and as such benefits from local geometric adaptiveness and broad convergence
properties. Due to its radial definition RVDE is continuous and computable in
linear time with respect to the dataset size. This amends for the main
shortcomings of previously studied VDEs, which are highly discontinuous and
computationally expensive. We provide a theoretical study of the modes of RVDE
as well as an empirical investigation of its performance on high-dimensional
data. Results show that RVDE outperforms other non-parametric density
estimators, including recently introduced VDEs.

##  1 INTRODUCTION

The problem of estimating a Probability Density Function (PDF) from a finite
set of samples lies at the heart of statistics and arises in several practical
scenarios ([10, 31]). Among density estimators, the non-parametric ones aim to
infer a PDF through a closed formula. Differently from parametric methods,
they do not require optimization and ideally provide an estimated PDF which is
simple, interpretable and computationally efficient. Two traditional examples
of non-parametric density estimators are the Kernel Density Estimator (KDE;
[16], [29]) and histograms ([14], [27]). KDE consists of a mixture of local
copies of a kernel around each datapoint while histograms partition the
ambient space into local cells (‘bins’) where the estimated PDF is constant.

Figure 1: An example of a density estimated via RVDE. The Voronoi tessellation
is depicted in solid gray. The estimated density is defined by the property
that its conical integral over the rays originating from the datapoints
(orange) is constant.

Both histograms and KDE suffer from bias due to the prior choice of a local
geometric structure i.e., the bins and the kernel respectively. This bias gets
exacerbated in high-dimensional ambient spaces. The reason is that datasets
grow exponentially in terms of geometric complexity, making a fixed simple
geometry unsuitable for estimating high-dimensional densities. This has led to
the introduction of the _Voronoi Density Estimator_ (VDE; [25]). VDE relies on
the geometric adaptiveness of Voronoi cells, which are convex polytopes
defined locally by the data ([24]). The PDF estimated by VDE is constant on
such cells, thus behaving as an adaptive version of histograms. Due to its
local geometric properties, VDE possesses convergence guarantees to the
ground-truth PDF which are more general than the ones of KDE.

The geometric benefits of VDE, however, come with a number of shortcomings.
First, the Voronoi cells and in particular their volumes are computationally
expensive to compute in high dimensions. Although this has been recently
attenuated by proposing Monte Carlo approximations ([28]), VDE falls behind
methods such as KDE in terms of computational complexity. Second, VDE
(together with its generalized version from [28] deemed CVDE) is highly
discontinuous on the boundaries of Voronoi cells. The estimated PDF
consequently suffers from large variance and instability with respect to the
dataset. This is again in contrast to KDE, which is continuous in its ambient
space.

In this work, we propose a novel non-parametric density estimator deemed
_Radial Voronoi Density Estimator_ (RVDE) which addresses the above
challenges. Similarly to VDE, RVDE integrates to a constant on Voronoi cells
and thus shares its local geometric advantages and convergence properties. In
contrast to VDE, RVDE is continuous and computable in linear time with respect
to the dataset size. The central idea behind RVDE is to define the PDF
radially from the datapoints so that the (conic) integral over the ray cast in
the corresponding Voronoi cell is constant (see Figure 1). This is achieved
via a ‘radial bandwidth’ which is defined implicitly by an integral equation.
Intuitively, the radial approach reduces the high-dimensional geometric
challenge of defining a Voronoi-based estimator to a one-dimensional problem.
This avoids the expensive volume computations of the original VDE and
guarantees continuity because of the fundamental properties of Voronoi
tessellations. Another important aspect of RVDE is its geometric distribution
of modes. We show that the modes either coincide with the datapoints or lie
along the edges of the Gabriel graph ([15]) depending on a hyperparameter
analogous to the bandwidth in KDE.

We compare RVDE with CVDE, KDE and the adaptive version of the latter in a
series of experiments. RVDE outperforms the baselines in terms of the quality
of the estimated density on a variety of datasets. Moreover, it runs
significantly faster and with lower sampling variance compared to CVDE. This
empirically confirms that the geometric and continuity properties of RVDE
translate into benefits for the estimated density in a computationally
efficient manner. We provide an implementation of RVDE (together with
baselines and experiments) in $C\texttt{++}$ at a publicly available
repository 111https://github.com/giovanni-marchetti/rvde. The code is
parallelized via the OpenCL framework and comes with a Python interface. In
summary our contributions include:

  * •

A novel density estimator (RVDE) based on the geometry of Voronoi
tessellations which is continuous and computationally efficient.

  * •

A complete study of the modes of RVDE and their geometric distribution.

  * •

An empirical investigation comparing RVDE to KDE (together with its adaptive
version) and previously studied VDEs.

##  2 RELATED WORK

###  2.1 Non-parametric Density Estimation

Non-parametric methods for density estimation trace back to the introduction
of histograms ([27]). Histograms have been extended by considering bin
geometries beyond the canonical rectangular one, for example triangular ([30])
and hexagonal ([4]) geometries. Another popular density estimator is KDE,
first discussed by [29] and [26]. The estimated density is a mixture of copies
of a priorly chosen distribution (‘kernel’) centered at the datapoints. KDE
has been extended to the multivariate case ([17, 7]) and has seen improvements
such as bandwidth selection methods ([21, 36]) and algorithms for adaptive
bandwidths ([37, 34]). Applications of KDE include estimation of traffic
incidents ([38]), of archaeological data ([1]) and of wind speed ([2]) to name
a few. As discussed in Section 1, both KDE and histograms suffer from lack of
geometric adaptiveness due to the choice of prior local geometries.

Another class of non-parametric methods are the orthogonal density estimators
([32, 22]). Those consist of choosing a discretized orthonormal basis of
functions and computing the coefficients of the ground-truth density via
Monte-Carlo integration over the dataset. When the basis is the Fourier one,
the estimator is referred to as ‘wavelet estimator’. The core drawback is that
orthogonal density estimators do not scale efficiently to higher dimensions.
When considering canonical tensor product bases the complexity grows
exponentially w.r.t. the dimensionality ([35]), making the estimator
unfeasible to compute.

###  2.2 Voronoi Density Estimators

The first Voronoi Density Estimator (VDE) has been pioneered by [25]. The
estimated density relies on Voronoi tessellations in order to achieve local
geometric adaptiveness. This is the main advantage of VDE over methods such as
KDE. The original VDE has seen applications to real-world densities such as
neurons in the brain ([12]), photons ([13]) and stars in a galaxy ([33]).
However, the method is not immediately extendable to high-dimensional spaces
because of unfeasible computational complexity of volumes and abundance of
unbounded Voronoi cells. This has been only recently amended by [28] by
introducing approximate numerical algorithms and by shaping of the density via
a kernel. In the present work, we aim to design an alternative version of the
original VDE which is continuous and does not rely on volume computations. The
resulting estimator is thus stable and computationally efficient while still
benefiting from the geometry of Voronoi tessellations.

##  3 BACKGROUND

In this section we recall the class of non-parametric density estimators which
we will be interested in throughout the present work. To this end, let
$P\subseteq\mathbb{R}^{n}$ be a finite set and consider the following central
notion from computational geometry.

######  Definition 1.

The _Voronoi cell_ 222Sometimes referred to as _Dirichlet cell_. of $p\in P$
is:

$C(p)=\\{x\in\mathbb{R}^{n}\ |\ \forall q\in P\ d(x,q)\geq d(x,p)\\}.$ (1)

The Voronoi cells are convex polytopes that intersect at the boundary and
cover the ambient space $\mathbb{R}^{n}$. The collection $\\{C(p)\\}_{p\in P}$
is referred to as _Voronoi tessellation_ generated by $P$. Note that although
the Voronoi tessellations are defined in an arbitrary metric space, the
resulting cells might be non-convex for distances different from the Euclidean
one. Since convexity will be crucial for the following constructions, we stick
to the Euclidean metric for the rest of the work.

We call _density estimator_ any mapping associating a probability density
function $f_{P}\in L^{1}(\mathbb{R}^{n})$ to a finite set
$P\subseteq\mathbb{R}^{n}$. The following class of density estimators
generalizes the original one by [25].

######  Definition 2.

A _Voronoi Density Estimator_ (VDE) is a density estimator $P\mapsto f_{P}$
such that for each $p\in P$:

$\int_{C(p)}f_{P}(x)\ \textnormal{d}x=\frac{1}{|P|}.$ (2)

VDEs stand out among density estimators for their geometric properties. This
is because the Voronoi cells are arbitrary polytopes that are adapted to the
local geometry of data. For VDEs all the Voronoi cells have the same estimated
probability, making such estimators locally adaptive from a geometric
perspective. This is reflected, for example, by the general convergence
properties of VDEs. The following is the main theoretical result from [28].

######  Theorem 3.1.

Let $P\mapsto f_{P}$ be a VDE and suppose that $P$ is sampled from a
probability density $\rho\in L^{1}(\mathbb{R}^{n})$ with support in the whole
$\mathbb{R}^{n}$. For $P$ of cardinality $m$ consider the probability measure
$\mathbb{P}_{m}=f_{P}\textnormal{d}x$ which is random in $P$. Then the
sequence $\mathbb{P}_{m}$ converges to $\mathbb{P}=\rho\textnormal{d}x$ in
distribution w.r.t. $x$ and in probability w.r.t. $P$. Namely, for any
measurable set $E\subseteq\mathbb{R}^{n}$ the sequence of random variables
$\mathbb{P}_{m}(E)$ converges in probability to the constant $\mathbb{P}(E)$.

In contrast, the convergence of other density estimators such as KDE requires
the kernel bandwidth to vanish asymptotically ([9]). The bandwidth vanishing
is necessary in order to amend for the local geometric bias inherent in KDE as
discussed in Section 1.

The following canonical construction of a VDE deemed Compactified Voronoi
Density Estimator (CVDE) is discussed by [28]. Given an integrable kernel $K:\
\mathbb{R}^{n}\times\mathbb{R}^{n}\rightarrow\mathbb{R}_{>0}$ the estimated
density is defined as

$f_{P}(x)=\frac{K(p,x)}{|P|\operatorname{\textnormal{Vol}}_{p}(C(p))}$ (3)

where $p$ is the closest point in $P$ to $x$ and
$\operatorname{\textnormal{Vol}}_{p}(C(p))=\int_{C(p)}K(p,y)\
\textnormal{d}y$. The latter volumes are approximated via Monte Carlo methods
since they become unfeasible to compute exactly as dimensions grow. The
resulting density inherits the same regularity as $K$ when restricted to each
Voronoi cell but jumps _discontinuously_ when crossing the boundary of Voronoi
cells (see Figure 3). Motivated by this, the goal of the present work is to
introduce a continuous and efficient VDE.

##  4 METHOD

###  4.1 Radial Voronoi Density Estimator

In this section we outline a general way of constructing a VDE with continuous
density function. Our central idea is to define the latter _radially_ w.r.t.
the datapoints $p\in P$. We start by rephrasing the integral over a Voronoi
cell (Equation 2) in spherical coordinates:

$\displaystyle\begin{split}&\int_{C(p)}f_{P}(x)\ \textnormal{d}x=\\\
=&\int_{\mathbb{S}^{n-1}}\underbrace{\int_{0}^{l(p+\sigma)}t^{n-1}f_{P}(p+t\sigma)\
\textnormal{d}t}_{\textnormal{Conical
Integral}}\textnormal{d}\sigma.\end{split}$ (4)

Here $\mathbb{S}^{n-1}\subseteq\mathbb{R}^{n}$ denotes the unit sphere and
$l(x)\in[0,+\infty]$ denotes the length of the segment contained in $C(p)$ of
the ray cast from $p$ passing through $x$ i.e.,

$\ l(x)=\sup\left\\{t\geq 0\ |\ p+t\frac{x-p}{d(x,p)}\in C(p)\right\\}.$ (5)

Figure 2: Illustration of the quantity $l(x)$ involved in the definition of
RVDE (Definition 3). The estimated density integrates to a constant over all
the infinitesimal cones (blue) in $C(p)$ originating from $p$.

We refer to Figure 2 for a visual illustration. Note that $l(x)$ is defined
for $x\not=p$ and is continuous in its domain since $l(x)=d(x,p)=d(x,q)$ for
$x\in C(p)\cap C(q)$.

We aim to solve Equation 4 by forcing the conical integral in Equation 4 to be
_constant_. To this end, we fix a continuous and strictly decreasing function
$K:\ \mathbb{R}_{>A}\rightarrow\mathbb{R}_{\geq 0}$ (a ‘kernel’) defined on a
half-line $\mathbb{R}_{>A}$, $A<0$, with the property that $t^{n-1}K(t)$ is
integrable on $\mathbb{R}_{>0}$. By an ansatz we look for a density in the
form

$f_{P}(x)=\frac{K\left(\beta(l(x))d(x,p)\right)}{\alpha|P|\textnormal{Vol}(\mathbb{S}^{n-1})}$
(6)

where $\alpha>0$ is a hyperparameter and $\beta:\
\mathbb{R}_{>0}\rightarrow\mathbb{R}$ is a function that we would like to
determine. The latter intuitively represents a _radial bandwidth_. The density
$f_{P}$ is continuous since the discontinuity of $l$ at $x=p$ is amended by
the vanishing of $d(x,p)$. Equation 2 is satisfied if for every $l>0$:

$\int_{0}^{l}t^{n-1}K(\beta(l)t)\ \textnormal{d}t=\alpha.$ (7)

Since $K$ is strictly decreasing, the above expression always has a unique
solution $\beta(l)>\frac{A}{l}$ assuming that $t^{n-1}K(t)$ is not integrable
around $A$. Such a guaranteed solution can be computed via any root-finding
algorithm and is continuous w.r.t. $l$. We provide an analysis of the function
$\beta$ and a discussion of the Newton-Raphson method for its computation in
Section 4.3.

The derivations above bring us to the following definition.

######  Definition 3.

Fix an $\alpha>0$ and a continuous function $K:\
\mathbb{R}_{>A}\rightarrow\mathbb{R}_{>0}$ with the domain bound $A<0$. Assume
the following:

  * •

$K(0)=1$,

  * •

$K$ is strictly decreasing,

  * •

$|t|^{n-1}K(t)$ is integrable around $+\infty$ but not integrable around $A$.

The _Radial Voronoi Density Estimator_ (RVDE) is the density estimator defined
by Equation 6 where $\beta$ is the function defined implicitly by Equation 7.

The following two standard families of kernels $K$ satisfy the above
requirements:

Exponential |  | Rational
---|---|---
$K(t)=e^{-t}$ |  | $K(t)=\frac{1}{(t+1)^{k}}$
(8)

where $k>n$. The domain bounds are $A=-\infty$ and $A=-1$ respectively. When
$n=1$ and $K$ is the exponential kernel, the function $\beta$ is closely
related to the Lambert $W$ function ([6]) via the expression:

$\beta(l)=\frac{1}{\alpha}+W\left(-\frac{l}{\alpha}e^{-\frac{l}{\alpha}}\right).$
(9)

We provide an empirical comparison between the two kernels from Equation 8 in
Section 5.3.

The intuition behind the hyperparameter $\alpha$ is that it controls the
trade-off between the amount of density concentrated around $P$ and away from
it (i.e., on the boundary of Voronoi cells). Indeed as $\alpha\to 0^{+}$ RVDE
tends (in distribution) to the discrete empirical measure over $P$ while as
$\alpha\to+\infty$ it tends to a measure concentrated on the boundary of
Voronoi cells. This can be deduced from Equation 7 since $\beta(l)$ tends to
$+\infty$ and to $A/l$ respectively and thus Equation 6 tends to $0$ for
$d(x,p)\not=0,l(x)$. This intuition around $\alpha$ will be corroborated by
Proposition 4.3, where we study how it controls the distribution of modes of
RVDE and consequently propose a heuristic selection procedure.

(a) KDE
([29])

(b) Original VDE
([25])

(c) CVDE
([28])

(d) RVDE
(Ours)

Figure 3: From left to right: heatmaps of KDE, of the two VDEs from the
literature and of our RVDE.

###  4.2 Computational Complexity and Sampling

We now discuss the computational cost of evaluating RVDE at a point
$x\in\mathbb{R}^{n}$. To begin with, the closest $p\in P$ to $x$ can be found
in logarithmic time w.r.t. $|P|$ by organizing $P$ in an efficient data
structure for nearest neighbor lookups such as a $k$-$d$ tree. Then $l(x)$ can
be computed in linear time via the following closed expression ([28]):

$l(x)=\min_{q\not=p,\ l^{q}(x)\geq 0}l^{q}(x)$ (10)

where

$l^{q}(x)=\frac{d(q,p)^{2}}{2\left\langle\frac{x-p}{d(x,p)},q-p\right\rangle}.$
(11)

The computational cost of evaluating $f_{P}(x)$ is thus linear w.r.t. $|P|$.
The remaining compute essentially reduces to solving Equation 7, which depends
on the integrator, the root-finder algorithm adopted and the desired
precision.

The formulation of RVDE enables a simple and efficient procedure for sampling
from the estimated density. In order to sample, one first chooses a $p\in P$
uniformly since $f_{P}$ integrates to $\frac{1}{|P|}$ on each Voronoi cell
(Equation 2). Since $t^{n-1}f_{P}(p+t\sigma)$ integrates to a constant on the
ray $r=\\{p+t\sigma\\}_{t\geq 0}\cap C(p)$ for every
$\sigma\in\mathbb{S}^{n-1}$, one then samples $\sigma$ uniformly from the
sphere. Finally one samples $t$ from the one-dimensional density $t^{n-1}K(t)$
restricted to the interval $[0,l(p+\sigma)]$. The computational complexity of
the latter step depends on the kernel as well as of the sampling method. The
result of the sampling is $p+t\sigma$. Because of the computational cost of
$l(p+\sigma)$, the sampling complexity of RVDE is linear w.r.t. $|P|$.

RVDE is more efficient than the VDE discussed by [28] (see the end of Section
3). The latter relies on Monte Carlo integration for numerical approximation
of volumes of Voronoi cells and has complexity $O(\Sigma|P|^{2})$ where
$\Sigma$ is the number of Monte Carlo samples. Compared to KDE, RVDE has the
same computational complexity (for both evaluation and sampling) while
retaining the geometric benefits of a VDE.

###  4.3 Study of $\beta$ and Modes

In this section we discuss qualitative properties and computational aspects of
the function $\beta$ defined implicitly by Equation 7 and consequently
characterize the modes of RVDE. We start by presenting an explicit expression
of the Newton-Raphson iteration for the computation of $\beta(l)$.

######  Proposition 4.1.

Fix $l>0$ and suppose $K\in C^{1}(\mathbb{R}_{>A})$ i.e., it is continuously
differentiable. Then the iteration $\beta_{m+1}$ of the Newton-Raphson method
for computing $\beta(l)$ by solving Equation 7 takes form:

$\displaystyle\beta_{m+1}=$ (12)
$=\beta_{m}+\frac{\beta_{m}}{n}\Biggl{(}1-\frac{l^{n}K(\beta_{m}l)-n\alpha}{l^{n}K(\beta_{m}l)-n\int_{0}^{l}t^{n-1}K(\beta_{m}t)\
\textnormal{d}t}\Biggl{)}.$

Moreover, if $K$ is convex then the Newton-Raphson method converges for any
initial value $\beta_{0}$ i.e., $\lim_{m\to\infty}\beta_{m}=\beta(l)$.

We refer to the Appendix for a proof. Note that the convexity assumption is
satisfied by both the kernels from Equation 8. Proposition 4.1 enables to
compute $\beta(l)$ and together with Section 4.2 provides all the algorithmic
details for implementing RVDE.

Next, we outline a qualitative study of the function $l\mapsto\beta(l)$.

######  Proposition 4.2.

The function $\beta:\ \mathbb{R}_{>0}\rightarrow\mathbb{R}$ is increasing, has
a zero at $l=(n\alpha)^{\frac{1}{n}}$ and has an horizontal asymptote:

$\lim_{l\to+\infty}\beta(l)=\left(\frac{1}{\alpha}\int_{0}^{\infty}t^{n-1}K(t)\
\textnormal{d}t\right)^{\frac{1}{n}}.$ (13)

Moreover if $K\in C^{1}(\mathbb{R}_{>A})$ then $\beta\in
C^{1}(\mathbb{R}_{>0})$ and it satisfies the differential equation:

$\left(l-\frac{n\alpha}{l^{n-1}K(\beta(l)l)}\right)\frac{\textnormal{d}\beta}{\textnormal{d}l}(l)=-\beta(l).$
(14)

We refer to the Appendix for a proof. As discussed in Section 4.1, $\beta$
generalizes the Lambert $W$ function. The properties and the differential
equation from Proposition 4.2 generalize their well-known instances for the
$W$ function ([6]).

We now focus on the study of modes. Our goal is to describe the modes of RVDE
completely. This is an advantage over density estimators such as KDE, where
the modes are challenging to describe and to compute approximately ([19, 5]).
Denote by $\varepsilon=(n\alpha)^{\frac{1}{n}}$ the zero of $\beta$.
Proposition 4.2 implies that for $x\in\mathbb{R}^{n}$, the density $f_{P}$
decreases radially w.r.t. $p$ in the direction of $x$ if $l(x)>\varepsilon$
and increases otherwise. This leads to the following result.

######  Proposition 4.3.

The modes of $f_{P}$ are classified as follows:

  1. (1)

$p\in P$ if $d(p,q)>2\varepsilon$ for every Voronoi cell $C(q)$ adjacent to
$C(p)$,

  2. (2)

$\frac{p+q}{2}$ for $p,q\in P$ if $\frac{p+q}{2}\in C(p)\cap C(q)$ and
$d(p,q)<2\varepsilon$,

  3. (3)

all points belonging to the segment $[p,q]$ for $p,q\in P$ if
$\frac{p+q}{2}\in C(p)\cap C(q)$ and $d(p,q)=2\varepsilon$.

We refer to the Appendix for a proof and to Figure 5 for an illustration.
Since $\varepsilon$ depends monotonically on the hyperparameter $\alpha$, the
latter controls the threshold for distances between adjacent points in $P$
below which the mode gets pushed away from such points towards the boundary of
the Voronoi cells. Intuitively, $\alpha$ determines the extent by which points
in $P$ are considered ‘isolated’ (i.e., a mode) or otherwise get ‘merged’ by
placing a mode between them.

(a) Gaussian

(b) Laplace

(c) Dirichlet

Figure 4: Comparison of the two kernels for RVDE (Equation 8) on three simple
distributions in $10$ dimensions.

An alternative geometric formulation of Proposition 4.3 is the following.
Consider the Gabriel graph of $P$ ([15]) containing an edge between $p$ and
$q$ iff $\frac{p+q}{2}\in C(p)\cap C(q)$ and discard all the edges of length
greater than $2\varepsilon$. The modes of RVDE are then associated with $(1)$
all isolated vertices, $(2)$ midpoints of edges and $(3)$ whole edges of
length $2\varepsilon$. Intuitively, the modes of RVDE are distributed
geometrically according to the truncated Gabriel graph.

Figure 5: Illustration of the modes of RVDE (red) together with the Gabriel
graph (black). Datapoints that not modes are marked in yellow.

This suggests a possible heuristic procedure for _hyperparameter selection_ of
$\alpha$. An option is to consider statistics of lengths of edges in the
Gabriel graph and choose $2\varepsilon$ (and thus $\alpha$) as a percentile.
The percentile we suggest is $\frac{|P|-1}{|E|}$ where $E$ denotes the set of
edges of the Gabriel graph. The intuition is that we wish to avoid modes
distributed in cycles. The number of cycles in the Gabriel graph is
$|P|-|E|+1$, from which our suggested percentile follows. This procedure
enables to select $\alpha$ automatically and we evaluate it empirically in
Section 5. However, it comes with a number of limitations. First, the
computational complexity of such a procedure is $O(|P|^{3})$ because of the
construction of the Gabriel graph, which is feasible but might become
expensive for large datasets. Another limitation is its independence from the
kernel $K$. The selection of $\alpha$ might be satisfying for some kernels but
not for others. In our empirical evaluation from Section 5 we show that for
the rational kernel the selected $\alpha$ is close to the optimal one in
practice, while for the exponential kernel the selection is further from
optimality.

##  5 EXPERIMENTS

Our empirical investigation is organized as follows. First we study RVDE on
its own by comparing the different choices of the kernel. We then compare RVDE
with other non-parametric density estimators on a variety of datasets.

###  5.1 Evaluation Metrics and Baselines

We evaluate all the density estimators $f_{P}$ via average log-likelihood on a
test set $P_{\textnormal{test}}$ i.e.,

$\frac{1}{|P_{\textnormal{test}}|}\sum_{x\in P_{\textnormal{test}}}\log
f_{P}(x).$ (15)

This measures whether the estimator assigns high density values to points
outside of $P$ but sampled from the same distribution. In order to empirically
evaluate the computational complexity, we additionally include runtimes for
all the experiments. Our implementations of all the considered density
estimators share the same programming framework and are parallelized to a
similar degree, making the raw runtimes a fair comparison. We perform
experiments on a machine with an AMD Ryzen 9 5950X 16-core CPU and a GeForce
RTX 3090 GPU.

We deploy the following non-parametric density estimators as baselines in the
experiments.

Kernel Density Estimator (KDE): given a (normalized) kernel
$K:\mathbb{R}^{n}\rightarrow\mathbb{R}_{\geq 0}$ the density is estimated as:

$f_{P}(x)=\frac{1}{|P|h^{n}}\sum_{p\in P}K\left(\frac{x-p}{h}\right)$ (16)

where $h$ is the bandwidth hyperparameter.

Adaptive Kernel Density Estimator (AdaKDE; [37]): a version of KDE where the
bandwidth $h_{p}$ depends on $p\in P$ and is smaller when data is denser
around $p$. Specifically, if $f_{P}(p)$ denotes the standard KDE estimate with
a global bandwidth $h$ then $h_{p}=h\lambda_{p}$, where:

$\lambda_{p}=(g/f_{P}(p))^{\frac{1}{2}},\quad\quad g=\prod_{q\in
P}f_{P}(q)^{\frac{1}{|P|}}.$ (17)

Compactified Voronoi Density Estimator (CVDE; [28]): the VDE described at the
end of Section 3. It depends on a kernel $K$ (together with a bandwidth) and
is discontinuous on the boundary of Voronoi cells.

###  5.2 Datasets

In our experiments we consider data of varying nature. This includes both
simple synthetic distributions and real-world datasets in high dimensions. For
the latter, we consider sound data ($n=21$) and image data ($n=100$). Our
datasets are the following.

Synthetic Datasets: datasets generated from a number of simple densities in
$n=10$ dimensions. Both $P$ and $P_{\textnormal{test}}$ contain $1000$ points
in all the cases. The densities we consider are: a standard Gaussian
distribution, a standard Laplace distribution, a Dirichlet distribution with
parameters $\alpha_{i}=\frac{1}{n+1}$ and a mixture of two Gaussians with
means $\mu_{1}=(-0.5,0,\cdots,0)$, $\mu_{2}=(0.5,0,\cdots,0)$ and standard
deviations $\sigma_{1}=0.1$, $\sigma_{2}=10$ respectively.

MNIST ([8]): a dataset consisting of $28\times 28$ grayscale images of
handwritten digits which are normalized in order to lie in $[0,1]^{28\times
28}$. In order to densify the data and obtain more meaningful estimates, we
downscale the images to resolution $10\times 10$. For each experimental run,
we sample half of the $60000$ training datapoints in order to evaluate the
variance of the estimation. The test set size is $10000$.

Anuran Calls ([11]): a dataset consisting of $7195$ calls from $10$ species of
frogs which are represented by $21$ normalized mel-frequency cepstral
coefficients in $[0,1]^{21}$. We retain $10\%$ of data for testing and sample
half of the training data at each experimental run.

(a) Single Gaussian

(b) Gaussian Mixture

(c) Anuran Calls

(d) MNIST

Figure 6: Comparison of the estimators as the bandwidth varies. All the
estimators implement the rational kernel.

###  5.3 Comparison of Kernels

Our first experiment consists of a comparison between the rational and the
exponential kernel for RVDE (Equation 8) on the synthetic datasets. In what
follows the exponent $k$ for the rational kernel is set to $k=n+1$ for
simplicity, where $n$ is the dimension of the ambient space of the dataset
considered (in this experiment, $n=10$).

The results are presented in Figure 4. The plot displays the test log-
likelihood (Equation 15) as the hyperparameter $\alpha$ varies. The latter is
scaled as $\alpha^{\frac{1}{n}}$ in order to be consistent with the
visualizations in the following section. The curves on the plot represent mean
and standard deviation (shaded areas) over $5$ experimental runs for $100$
bandwidths. The additional vertical lines correspond to the value of the
hyperparameter selection heuristic discussed at the end of Section 4.3. As can
be seen, the performance of the rational kernel is more stable w.r.t. the
hyperparameter $\alpha$. The exponential kernel, however, achieves a slightly
higher test score with its best hyperparameter on the Gaussian and Laplace
datasets. Note that the heuristically chosen $\alpha$ aligns well with the one
that achieves the best performance for the rational kernel, but is misaligned
for the exponential one. We conclude that the rational kernel is generally a
better option unless an extensive hyperparameter search is performed. In what
follows we consequently stick to the rational kernel for RVDE.

###  5.4 Comparison with Baselines

In our main experiment we compare the performance of RVDE with the baselines
described in Section 5.1. We consider the test log-likelihood (Equation 5.1),
the standard deviation of the latter and the runtimes. In order to make the
comparison as fair as possible, all the estimators are implemented with the
rational kernel. We found out that the performances drop with the more
standard Gaussian kernel (which does not apply to RVDE). We include the
results with both the Gaussian kernel and the exponential kernel in the
Appendix.

Table 1: Average runtimes (in seconds) per one full train-test run with fixed
bandwidth. RVDE is highlighted in blue.

| RVDE | CVDE | KDE | AdaKDE
---|---|---|---|---
Gaussian | 0.0376 | 0.265 | 0.0340 | 0.266
Anuran Calls | 0.0581 | 0.490 | 0.0787 | 0.870
MNIST | 17.4 | 408 | 12.5 | 75.0

The plot in Figure 6 displays the (test) log-likelihood for all the estimators
as the bandwidth hyperparameter $h$ varies (see the definition of the
baselines in Section 5.1). In order to compare RVDE on the same scale as the
other estimators, we convert $h$ to $\alpha$ via:

$\alpha=\int_{0}^{\infty}K\left(\frac{t}{h}\right)\
\textnormal{d}t=h^{n}\int_{0}^{\infty}K(t)\ \textnormal{d}t.$ (18)

As can be seen, RVDE outperforms the baselines (each with its respective best
bandwidth) in all the cases considered. The margin between RVDE and the
baselines is especially evident on the more complex and high-dimensional
datasets (Anuran Calls and MNIST). This confirms that the geometric benefits
and the continuity properties of RVDE translate into better estimates for
densities of different nature and increasing dimensionality.

Table 1 reports the average runtime for an experimental run (with a single
fixed bandwidth) for each estimator. RVDE outperforms the CVDE as well as
AdaKDE by an extremely large margin. KDE achieves comparable runtimes to RVDE:
it is slightly faster on Gaussian and MNIST while it is slightly slower on
Anuran Calls. This confirms empirically the discussion from Section 4.2: RVDE
is significantly more efficient than CVDE and has the same (asymptotic)
complexity as KDE.

Table 2: Standard deviations of the (test) log-likelihood over $5$
experimental runs. RVDE is highlighted in blue. Each estimator is considered
with its best bandwidth.

| RVDE | CVDE | KDE | AdaKDE
---|---|---|---|---
Gaussian | 0.788 | 0.843 | 0.572 | 0.572
Anuran Calls | 1.170 | 1.253 | 1.152 | 1.152
MNIST | 5.507 | 5.767 | 5.735 | 5.735

Table 2 separately reports the standard deviation of the log-likelihood
(averaged over $P_{\textnormal{test}}$) w.r.t. the dataset sampling. For each
estimator, we consider its best bandwidth according to the results from Figure
6. We first observe that RVDE achieves lower standard deviation than CVDE on
all the datasets. This corroborates the hypothesis that the continuity of RVDE
results in more stable estimates than those obtained by the highly-
discontinuous CVDE. KDE and AdaKDE achieve the lowest standard deviations on
the Gaussian and Anuran Calls datasets. This is likely due to the smoothness
of such estimators and again confirms the benefit of regularity biases in
terms of stability. However, on the most complex dataset considered (MNIST)
RVDE outperforms the baselines. This suggests that for articulated densities
the biases of geometric nature become more beneficial than generic biases such
as smoothness.

##  6 CONCLUSIONS AND FUTURE WORK

In this work we introduced a non-parametric density estimator (RVDE)
benefiting from the geometric properties of Voronoi tessellations while being
continuous and computationally efficient. We provided both theoretical and
empirical investigations of RVDE.

An interesting line for future investigation is to explore the radial
construction of RVDE on Riemannian manifolds beyond the Euclidean space. In
this generality the rays correspond to geodesics defined via the exponential
map of the given manifold. A variety of Riemannian manifolds arise in
statistics and machine learning. For example, data on spheres are the object
of study of directional statistics ([20]), hyperbolic spaces are routinely
deployed to represent hierarchical data ([23]) and complex projective spaces
correspond to Kendall shape spaces from computer vision ([18]). Those areas of
research can potentially benefit from the geometric characteristics and the
computational efficiency of an extension of RVDE to Riemannian manifolds.

##  7 Acknowledgements

This work was supported by the Swedish Research Council, Knut and Alice
Wallenberg Foundation and the European Research Council (ERC-BIRD-884807).

## References

  * [1] Michael J. Baxter, Christian C. Beardah and Richard V.. Wright  “Some archaeological applications of kernel density estimates”  In _Journal of Archaeological Science_ 24.4 Elsevier, 1997, pp. 347–354
  * [2] Hu Bo, Li Yudun, Yang Hejun and Wang He  “Wind speed model based on kernel density estimation and its application in reliability assessment of generating systems”  In _Journal of Modern Power Systems and Clean Energy_ 5.2 Springer, 2017, pp. 220–227
  * [3] Stephen Boyd, Stephen P. Boyd and Lieven Vandenberghe  “Convex optimization”  Cambridge university press, 2004
  * [4] Daniel B. Carr, Anthony R. Olsen and Denis White  “Hexagon mosaic maps for display of univariate and bivariate geographical data”  In _Cartography and Geographic Information Systems_ 19.4 Taylor & Francis, 1992, pp. 228–236
  * [5] Dorin Comaniciu and Peter Meer  “Mean shift: A robust approach toward feature space analysis”  In _IEEE Transactions on pattern analysis and machine intelligence_ 24.5 IEEE, 2002, pp. 603–619
  * [6] Robert M. Corless et al.  “On the LambertW function”  In _Advances in Computational mathematics_ 5.1 Springer, 1996, pp. 329–359
  * [7] Khosrow Dehnad  “Density estimation for statistics and data analysis”  Taylor & Francis, 1987
  * [8] Li Deng  “The mnist database of handwritten digit images for machine learning research”  In _IEEE Signal Processing Magazine_ 29.6 IEEE, 2012, pp. 141–142
  * [9] L.. Devroye and T.. Wagner  “The L1 convergence of kernel density estimates”  In _The Annals of Statistics_ JSTOR, 1979, pp. 1136–1139
  * [10] Peter J. Diggle  “Statistical analysis of spatial and spatio-temporal point patterns”  CRC press, 2013
  * [11] Dheeru Dua and Casey Graff  “UCI Machine Learning Repository”, 2017  URL: http://archive.ics.uci.edu/ml
  * [12] Charles Duyckaerts, Gilles Godefroy and Jean-Jacques Hauw  “Evaluation of neuronal numerical density by Dirichlet tessellation”  In _Journal of neuroscience methods_ 51.1 Elsevier, 1994, pp. 47–69
  * [13] H. Ebeling and G. Wiedenmann  “Detecting structure in two dimensions combining Voronoi tessellation and percolation”  In _Physical Review E_ 47.1 APS, 1993, pp. 704
  * [14] David Freedman and Persi Diaconis  “On the histogram as a density estimator: L 2 theory”  In _Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete_ 57.4 Springer, 1981, pp. 453–476
  * [15] K. Gabriel and Robert R. Sokal  “A new statistical approach to geographic variation analysis”  In _Systematic zoology_ 18.3 Society of Systematic Zoology, 1969, pp. 259–278
  * [16] Artur Gramacki  “Nonparametric kernel density estimation and its computational aspects”  Springer, 2018
  * [17] Alan Julian Izenman  “Review papers: Recent developments in nonparametric density estimation”  In _Journal of the american statistical association_ 86.413 Taylor & Francis, 1991, pp. 205–224
  * [18] Christian Peter Klingenberg  “Walking on Kendall’s shape space: understanding shape spaces and their coordinate systems”  In _Evolutionary Biology_ 47.4 Springer, 2020, pp. 334–352
  * [19] Jasper C.. Lee et al.  “Finding an approximate mode of a kernel density estimate”  In _29th Annual European Symposium on Algorithms (ESA 2021)_ , 2021  Schloss Dagstuhl-Leibniz-Zentrum für Informatik
  * [20] Kanti V. Mardia and Peter E. Jupp  “Directional statistics”  John Wiley & Sons, 2009
  * [21] J.. Marron  “A comparison of cross-validation techniques in density estimation”  In _The Annals of Statistics_ JSTOR, 1987, pp. 152–162
  * [22] Elias Masry  “Multivariate probability density estimation by wavelet methods: Strong consistency and rates for stationary time series”  In _Stochastic processes and their applications_ 67.2 Elsevier, 1997, pp. 177–193
  * [23] Maximillian Nickel and Douwe Kiela  “Poincaré embeddings for learning hierarchical representations”  In _Advances in neural information processing systems_ 30, 2017
  * [24] Atsuyuki Okabe, Barry Boots, Kokichi Sugihara and Sung Nok Chiu  “Spatial tessellations: concepts and applications of Voronoi diagrams”  John Wiley & Sons, 2009
  * [25] J.. Ord  “How many trees in a forest”  In _Mathematical Scientist_ 3, 1978, pp. 23–33
  * [26] Emanuel Parzen  “On estimation of a probability density function and mode”  In _The annals of mathematical statistics_ 33.3 JSTOR, 1962, pp. 1065–1076
  * [27] Karl Pearson  “Contributions to the mathematical theory of evolution”  In _Philosophical Transactions of the Royal Society of London. A_ 185 JSTOR, 1894, pp. 71–110
  * [28] Vladislav Polianskii et al.  “Voronoi density estimator for high-dimensional data: Computation, compactification and convergence”  In _Uncertainty in Artificial Intelligence_ , 2022, pp. 1644–1653  PMLR
  * [29] Murray Rosenblatt  “Remarks on Some Nonparametric Estimates of a Density Function”  In _The Annals of Mathematical Statistics_ 27.3 Institute of Mathematical Statistics, 1956, pp. 832–837  URL: http://www.jstor.org/stable/2237390
  * [30] David W. Scott  “A note on choice of bivariate histogram bin shape”  In _Journal of Official Statistics_ 4.1 Statistics Sweden (SCB), 1988, pp. 47
  * [31] David W. Scott  “Multivariate density estimation: theory, practice, and visualization”  John Wiley & Sons, 2015
  * [32] Marina Vannucci  “Nonparametric density estimation using wavelets”  Institute of Statistics & Decision Sciences, Duke University, 1995
  * [33] Iryna Vavilova, Andrii Elyiv, Daria Dobrycheva and Olga Melnyk  “The Voronoi tessellation method in astronomy”  In _Intelligent Astrophysics_ Springer, 2021, pp. 57–79
  * [34] Christiaan Maarten Walt and Etienne Barnard  “Variable kernel density estimation in high-dimensional feature spaces”  In _Thirty-first AAAI conference on artificial intelligence_ , 2017
  * [35] G.. Walter  “Estimation with wavelets and the curse of dimensionality”  In _Manuscript, Department of Mathematical Sciences, University of Wisconsin-Milwaukee_ , 1995
  * [36] Matt P. Wand and M. Jones  “Multivariate plug-in bandwidth selection”  In _Computational Statistics_ 9.2 Heidelberg: Physica-Verlag,[1992-, 1994, pp. 97–116
  * [37] Bin Wang and Xiaofeng Wang  “Bandwidth selection for weighted kernel density estimation”  In _arXiv preprint arXiv:0709.1616_ , 2007
  * [38] Zhixiao Xie and Jun Yan  “Kernel density estimation of traffic accidents in a network space”  In _Computers, environment and urban systems_ 32.5 Elsevier, 2008, pp. 396–406

An Efficient and Continuous Voronoi Density Estimator:
Supplementary Materials

## PROOFS OF RESULTS FROM SECTION 4.3

######  Proposition 7.1.

Fix $l>0$ and suppose $K\in C^{1}(\mathbb{R}_{>A})$. Then the iteration
$\beta_{m+1}$ of the Newton-Raphson method for computing $\beta(l)$ by solving
Equation 7 takes form:

$\beta_{m+1}=\beta_{m}\left(1+\frac{1}{n}\left(1-\frac{l^{n}K(\beta_{m}l)-n\alpha}{l^{n}K(\beta_{m}l)-n\int_{0}^{l}t^{n-1}K(\beta_{m}t)\
\textnormal{d}t}\right)\right).$ (19)

Moreover, if $K$ is convex then the Newton-Raphson method converges for any
initial value $\beta_{0}$ i.e., $\lim_{m\to\infty}\beta_{m}=\beta(l)$.

###### Proof.

Consider

$F(\beta)=\int_{0}^{l}t^{n-1}K(\beta t)\ \textnormal{d}t-\alpha.$ (20)

The iteration of the Newton-Rhapson method for solving $F(\beta)=0$ takes
form:

$\beta_{m+1}=\beta_{m}-\frac{F(\beta_{m})}{\frac{\textnormal{d}F}{\textnormal{d}\beta}(\beta_{m})}.$
(21)

Via integration by parts we obtain:

$\frac{\textnormal{d}F}{\textnormal{d}\beta}(\beta)=\int_{0}^{l}t^{n}\frac{\textnormal{d}K}{\textnormal{d}t}(\beta
t)\ \textnormal{d}t=\frac{1}{\beta}\left(l^{n}K(\beta
l)-n\int_{0}^{l}t^{n-1}K(\beta t)\ \textnormal{d}t\right).$ (22)

Equation 19 follows then from Equation 21 by elementary algebraic
manipulations. The convergence guarantee follows from the fact that if $K$ is
convex then $F$ is easily seen to be convex as well. The Newton-Raphson method
is well-known to be convergent for convex functions ([3]). ∎

######  Proposition 7.2.

The function $\beta:\ \mathbb{R}_{>0}\rightarrow\mathbb{R}$ is increasing, has
a zero at $l=(n\alpha)^{\frac{1}{n}}$ and has an horizontal asymptote:

$\lim_{l\to+\infty}\beta(l)=\left(\frac{1}{\alpha}\int_{0}^{\infty}t^{n-1}K(t)\
\textnormal{d}t\right)^{\frac{1}{n}}.$ (23)

Moreover if $K\in C^{1}(\mathbb{R}_{>A})$ then $\beta\in
C^{1}(\mathbb{R}_{>0})$ and it satisfies the differential equation:

$\left(l-\frac{n\alpha}{l^{n-1}K(\beta(l)l)}\right)\frac{\textnormal{d}\beta}{\textnormal{d}l}(l)=-\beta(l).$
(24)

###### Proof.

The claim on the monotonicity of $\beta$ follows directly from its definition
(Equation 7) and the hypothesis that $K$ is decreasing. In order to compute
its zero, note that $\beta(l)=0$ implies
$\alpha=\int_{0}^{l}K(0)t^{n-1}=\frac{l^{n}}{n}$ and thus
$l=(n\alpha)^{\frac{1}{n}}$. For the asymptote note that for $l=+\infty$
Equation 7 becomes by a change of variables:

$\int_{0}^{\infty}t^{n-1}K(\beta(+\infty)t)\
\textnormal{d}t=\frac{1}{\beta(+\infty)^{n}}\int_{0}^{\infty}t^{n-1}K(t)\
\textnormal{d}t=\alpha.$ (25)

Lastly, in order to obtain the differential equation for $\beta$ we
differentiate Equation 7 on both sides and get:

$\begin{split}0=\frac{\textnormal{d}}{\textnormal{d}l}\int_{0}^{l}t^{n-1}K(\beta(l)t)\
\textnormal{d}t&=l^{n-1}K(\beta(l)l)+\int_{0}^{l}t^{n-1}\frac{\textnormal{d}}{\textnormal{d}l}K(\beta(l)t)\
\textnormal{d}t\\\
&=l^{n-1}K(\beta(l)l)+\frac{\textnormal{d}\beta}{\textnormal{d}l}(l)\int_{0}^{l}t^{n}\frac{\textnormal{d}K}{\textnormal{d}t}(\beta(l)t)\
\textnormal{d}t\\\
&=l^{n-1}K(\beta(l)l)+\frac{\textnormal{d}\beta}{\textnormal{d}l}(l)\frac{l^{n}K(\beta(l)l)-n\alpha}{\beta(l)}\end{split}$
(26)

where in the first identity we deployed the (distributional) Leibniz rule
while in the last one we deployed integration by parts.

∎

######  Proposition 7.3.

The modes of $f_{P}$ are as follows:

  1. (1)

$p\in P$ if $d(p,q)>2\varepsilon$ for every Voronoi cell $C(q)$ adjacent to
$C(p)$,

  2. (2)

$\frac{p+q}{2}$ for $p,q\in P$ if $\frac{p+q}{2}\in C(p)\cap C(q)$ and
$d(p,q)<2\varepsilon$,

  3. (3)

all points belonging to the segment $[p,q]$ for $p,q\in P$ if
$\frac{p+q}{2}\in C(p)\cap C(q)$ and $d(p,q)=2\varepsilon$.

###### Proof.

Pick $p\in P$. If $p$ satisfies the hypothesis of the first claim then
$l(x)>\varepsilon$ for every $x\in C(p)$ and thus $\beta(l(x))>0$ by
Proposition 7.2. Since $K$ is decreasing, $f_{P}$ decreases radially w.r.t.
$p$ in $C(p)$ and the first claim follows. If $p$ does not satisfy the
hypothesis of the first claim then $\beta(l(x))\leq 0$ for some $x\in C(p)$.
With the exception of the case $\beta(l)=0$, the modes lie then on the
boundary and are of the form $K(\beta(l)l)$ up to a multiplicative constant.
The function $\beta(l)l$ is increasing in $l$ since by appealing to
Proposition 7.2 we can compute its derivative:

$\frac{\textnormal{d}\beta(l)l}{\textnormal{d}l}=\beta(l)+l\frac{\textnormal{d}\beta(l)}{\textnormal{d}l}=\beta(l)\frac{n\alpha}{n\alpha-l^{n}K(\beta(l)l)}\geq
0.$ (27)

Since $l(x)$ has local minima at midpoints of segments connecting points in
$P$, $K(\beta(l)l)$ is locally maximized therein and the second claim follows.
In the hypothesis of the third claim $\beta$ vanishes on the segment and the
density is thus constant.

∎

## ADDITIONAL EXPERIMENTS

In this section we report additional experimental results complementing the
ones in the main of the manuscript. For completeness, we evaluate the density
estimators on different kernels. Figure 8 displays comparative results for all
the estimators with the exponential and the Gaussian kernel (note that the
latter does not apply to RVDE). Moreover, we experiment with different
dimensions and evaluation metrics other than average log-likelihood. This is
possible only on a synthetic dataset where the dimension can vary and where
the ground-truth density $\rho$ is known. The latter is necessary for the
metric considered. Figure 8 displays a comparison on a high-dimensional
Gaussian mixture ($n=30$) as well as a comparison on the Gaussian mixture as
in Section 5 ($n=10$) where the evaluation metric is the empirical _Hellinger
distance_ on the test set:

$\frac{1}{2|P_{\rm test}|}\sum_{x\in P_{\rm
test}}\left(f_{P}(x)^{\frac{1}{2}}-{\rho}(x)^{\frac{1}{2}}\right)^{2}.$ (28)

(a) Single Gaussian

(b) Gaussian Mixture

(c) Anuran Calls

(d) MNIST

Exponential

Gaussian

Figure 7: Comparison of the estimators with the exponential and Gaussian
kernel as the bandwidth varies.

(e) $n=30$

(f) Hellinger Distance

Figure 8: Comparison of the estimators on a $30$-dimensional Gaussian mixture
(left) and on a $10$-dimensional Gaussian mixture with the Hellinger distance
as a metric (right).