Voronoi Boundary Classification: A High-Dimensional Geometric Approach
via Weighted Monte Carlo Integration

Vladislav Polianskii 1 Florian T. Pokorny 1

Abstract
Voronoi cell decompositions provide a classicalavenue to classification. Typical approaches however only utilize point-wise cell-membership information by means of nearest neighbor queries
and do not utilize further geometric information
about Voronoi cells since the computation ofVoronoi diagrams is prohibitively expensive inhigh dimensions. We propose a Monte-Carlo integration based approach that instead computes aweighted integral over the boundaries of Voronoi
cells, thus incorporating additional information
about the Voronoi cell structure. We demonstrate
the scalability of our approach in up to 3072 dimensional spaces and analyze convergence basedon the number of Monte Carlo samples and choiceof weight functions. Experiments comparing ourapproach to Nearest Neighbors, SVM and Random Forests indicate that while our approach
performs similarly to Random Forests for large
data sizes, the algorithm exhibits non-trivial datadependent performance characteristics for smaller
datasets and can be analyzed in terms of a geometric confidence measure, thus adding to therepertoire of geometric approaches to classification while having the benefit of not requiring any
model changes or retraining as new training samples or classes are added.

1. Introduction and Related Work
The problem of classifying a set of n points {x1 , . . . , xn } ind-dimensional Euclidean space, each labelled with one of kclasses has been tackled with a large variety of approaches
in recent decades. In many application fields, Deep Neuralnetworks are today leading the field in terms of classifica1

Division of Robotics, Perception and Learning, EECS, KTH
Royal Institute of Technology, Stockholm, Sweden. Correspondence to: Vladislav Polianskii <vpol@kth.se>.
Proceedings of the 36 th International Conference on Machine
Learning, Long Beach, California, PMLR 97, 2019. Copyright2019 by the author(s).

tion rates, particularly in the area of image classification(Krizhevsky et al., 2012). A draw-back of these approaches
to date is the difficulty of reasoning about the learned representations in terms of convergence guarantees, stability
and the complex structure of learned embeddings leading tochallenging phenomena such as adverserial examples (Kurakin et al., 2016). In contrast, a benefit of purely geometricapproaches such as nearest neighbors – while not always
having optimal performance – is their simplicity and asymptotic convergence under suitable assumptions (Döring et al.,2017). A further advantage is the lack of dependence on initialization since no stochastic gradient descent or complex
variational inference procedure is required. This work contributes to this category of optimization-free, non-Bayesianand purely geometric approaches to classification with afocus on simpliticy and mathematical interpretability.
Recent works such as (Sharif Razavian et al., 2014) havedemonstrated that classical methods such as SVM classifiers
can provide competitive results when combined with representations determined by a neural network. Similarly, metriclearning techniques (Weinberger & Saul, 2009) may be used
to first determine an appropriate metric and/or embedding
before applying geometric methods such as a nearest neighbor classifier. The Voronoi cell based approach proposed
here can hence also be viewed as a component in the toolbox
that may be combined with Bayesian or Neural Network
based methods. In the present paper, we focus on the definition, algorithms, mathematical properties and an initialstand-alone evaluation of our approach.
While nearest-neighbor classifiers (Biau & Devroye, 2015)
make a classification decision by testing which training data
Voronoi cell a test data point t ∈ Rd belongs to, this membership test can be performed efficiently with approximate
geometric lookup datastructures (Kushilevitz et al., 2000)
and does not require a computation of the full geometryof the Voronoi cells themselves. Our approach is instead
focused on the Voronoi cell V (t) around the test data pointt and balances the influence of neighboring training data on
the classification decision by means of a weighted integral
over the Voronoi boundary. A related approach consistsof the Natural Neighbor regression method (Sibson, 1981),
where volumes of intersections of Voronoi cells before andVoronoi Boundary Classification: A High-Dimensional Geometric Approach via Weighted Monte Carlo Integration

after insertion of a test-point are used as weights for interpolation. By considering a ranking of these weights thismethod could be considered as a related classification approach, but since current implementations rely on numericalintegration and an explicit computation of the Voronoi cellstructure, these are not applicable for high-dimensional datain their current form.
The key contributions of this work are: 1) We introduce
and provide a convergence analysis of a conceptually simple Voronoi-based classifier based on a weighted integral
formulation over boundary faces of V (t). 2) We show how
to overcome the computational complexity challenge ofworking with Voronoi cells in high dimensions by avoiding explicit calculation of the Voronoi structure. This isachieved by means of a Monte Carlo integral approximation
method and, to the best of our knowledge, is the first time
an integral over Voronoi cell geometry has been proposed
and used for classification in high dimensions. 3) We evaluate our approach relative to Nearest Neighbors, SVM and
Random Forests on 3 datasets, including CIFAR-10 and
MNIST.

2. Methodology and Algorithm
We consider a test data point t ∈ Rd and a labeled training
dataset X = {x1 , . . . , xn } ⊂ Rd where each data pointxi ∈ Rd lies in Euclidean space and each xi is associated toa unique class ci ∈ {1, . . . , k}. We assume furthermore thatthe points in X are in general position, which for data withuniform noise happens with probability one and which canalternatively be enforced by an arbitrarily small perturbationof the input data. For any x ∈ Rd , we denote the Voronoi
cell at x with respect to {t, x1 , . . . , xn } by V (x) = {y ∈Rd : ky − xk ≤ ky − xj k for all i ∈ {1, . . . , n} and ky −xk ≤ ky − tk}.

We seek to determine a classification label for t and consider
the Voronoi cell V (t) as the basic geometric object uponwhich a classification decision will be made. Here, we thinkof each boundary face of V (t) as being colored by a uniquecolor corresponding to the class label of the neighboring xi(see Fig.1). In intuitive terms, our basic motivation for the
proposed method is that an observer placed at t may decideupon the classification label for t by ranking the relative influence of the observed colored faces of V (t) by weighting
the distances between t and points on the boundary and thetotal area of the observed colored faces. I.e. the class of t isdetermined by “how much of a given color is observed fromt, weighted by distance”. We introduce a weight function
w : Rd → R+ that decays with distance from t and our
proposed method utilizes a w-weighted integral over theboundary faces of V (t) to achieve the above intuitive goal.
Crucially, we will show in the following that this formulation can be implemented without requiring the computation

Figure 1. Planar classification problem with two classes of training
data (purple / orange) as well as test data point t in black. Our
approach ranks the influence of each neighboring class on t, byapproximating a weighted integral over the Voronoi cell boundaryof t for each color.

of the full Voronoi cell structure, which is of complexitydO(nd 2 e ) (Aurenhammer, 1991) and thus infeasible in highdimensions.
Definition 1. Consider a labeled training dataset{(x1 , c1 ), . . . , (xn , cn )}, where xi ∈ Rd and ci ∈{1, . . . , k} for all i ∈ {1, . . . n}. For any integrable function w : Rd → R+ , t ∈ Rd and i ∈ {1, . . . , k}, the(weighted) Voronoi Boundary Rank of class i at t is denotedby ri (t) and given byri (t) =X Z
j:cj =iw(x)dVol,

V (t)∩V (xj )where dVol denotes the Volume form on V (t) ∩ V (xj )induced by the Euclidean metric in Rd . The (weighted)
Voronoi Boundary Classifier is given by V BC(t) =argmaxi∈{1,...,k} ri (t).
Observe that, in the above, each face V (t) ∩ V (xj ) is d − 1
dimensional (or empty) and may be unbounded. As a result,for w to be integrable over this domain, we consider w thatdecay with distance to t. As we shall explore later, exponentially decaying weight functions constitute a suitable choice,for example. We now turn to the problem of evaluating the
integral ri (t). We approach this by first changing variables
and performing the integration over a unit sphere instead.
Computation of Voronoi Boundary Rank
Lemma 1 (Change of Variables and Monte Carlo). Let
P = V (t) ∩ V (xj ). If P 6= ∅, P lies in some hyperplaneΠ = {x ∈ Rd : hx, ni = c}, where n is of unit norm.
Denote by Sd−1 the d − 1 dimensional unit sphere centeredat the origin and denote by fj the map which maps m ∈Sd−1 to to the intersection between the ray starting at t
in direction m with P if this intersection exists. Let Ij =

Voronoi Boundary Classification: A High-Dimensional Geometric Approach via Weighted Monte Carlo IntegrationRV (t)∩V (xj )w(x)dVol, thenZ
Ij =fj−1 (V (t)∩V (xj ))w(fj (m))

2δ
kfj (m) − tkd−1
dVolSd−1|hm, ni|(1)

where dVolSd−1 denotes the standard induced volume formof the sphere in Euclidean space. Secondly, consider a sequence {mi }T1 of uniform samples on fj−1 (V (t)∩V (xj )) ⊂Sd−1 . We haveT1X
kfj (mi ) − tkd−1
w(fj (mi ))T →∞ T|hmi , ni|
i=1

Ij = Vj lim(2)where Vj denotes the volume of the set fj−1 (V (t) ∩ V (xj ))with respect to the standard induced metric on the unit
sphere in Euclidean space.
Proof. The first part above follows from the general integral change of variable formula and a consideration of
Sylvester’s determinant theorem to simplify the resulting
expression, while the second part follows from the standard
Monte Carlo integration convergence (Davis & Rabinowitz,2007). See the supplementary material for additional details.
We note that the above sampling based approximation of
Ij still requires us to 1) perform uniform sampling on
fj−1 (V (t) ∩ V (xj )) ⊂ Sd−1 and 2) to determine its volume Vj . We note that this can be achieved by means ofrejection sampling by noting the following lemma.
dd−1

Lemma 2. Let t ∈ R and m ∈ S

and definekx − tk2
2 hm, x − ti∗∗min {lm(xi ) | lm(xi ) > 0} .

∗
lm
(x) =

face(m) = argi∈{1,...,n}∗
If there is no i ∈ {1, . . . , n} such that lm(xi ) > 0, then the
ray R starting at t and in direction m is fully contained inV (t) (and V (t) is unbounded in direction m starting at t).
Otherwise m lies in fj−1 (V (t) ∩ V (xface(m) )).

Proof. See the supplementary material for a proof.
Given Lemma 2, we can now uniformly sample T
samples ST = {m1 , . . . , mT } on Sd−1 and count the
number Tj (ST ) = #{m ∈ ST : face(m) = j} ofsamples mapped to V (t) ∩ V (xj ). This yields Vj =T (S )I
Vol(Sn−1 ) limT →∞ j T T . We now have: Vol(Sjn−1 ) =PT
kf (m )−tkd−1
limT →∞ T1 i=1 1face(mi )=j w(f (mi )) j |hmi i ,ni|
where 1face(mi )=j denotes the indicator function thatis 1 if face(mi ) = j and zero otherwise. Since

δ
0
−δ
−2δ

−2ε −ε0δ2δFigure 2. A counter-example for a bounded weight function.

Pri (t) = j:cj =i Ij , this provides the basis for computing
the Voronoi Classification scheme we propose here.
Observe that we in particular did not require an explicit
computation of the boundary face geometry or volumesof V (t). Additionally, a notable quality of the algorithm
is that the Monte Carlo integration rate of convergence isO( √1T ) in arbitrary dimension (Doucet et al., 2001).
Weight functions We now come to a discussion of suitable weight functions w for our proposed method. Observe
that, since Voronoi Cells may be unbounded, we may encounter an integral over an unbounded domain in the definition of the Voronoi rank. To ensure that all integrals
are finite we hence require that w decays with distance tothe test data point t. In the following experiments, we infact only consider w(x) = ŵ(kx − tk) for some integrable
function ŵ : R+ → R+ that is monotonically decreasing.

An obvious choice may be a weight function of the form
−2 2ŵ(z) = e−0.5σ z which does indeed result in finite integrals, but does - as we shall now explain - have the property
of favoring volume over distance in the relative contributionof the boundary of V (t). To illustrate this, consider points in
R2 , placed as shown on Figure 2: a test point t is placed at
the origin, one point of the first class at coordinates (−2ε, 0)
for some positive ε and 3 points of the second class at coordinates (0, 2δ), (0, −2δ) and (2δ, 0) for some positive δ.
In this√ setting, it is easy to show that r1 (t) ≤ ŵ(ε) · 2δ andŵ(δ 2)(4δ + 2ε) ≤ r2 (t). If we decrease ε towards 0, thepoint of the first class starts to move closer to the test point.
However, a simple analysis shows that if δ < √12 ŵ−1 (0.5),
then no matter how close the first point is moved towards
the origin, r2 (t) > r1 (t) and the classification chooses thesecond (orange) class - in effect the larger volume of theorange boundary component of V (t) always dominates theclassification.
This behavior may be considered unnatural if a balance between distance and volume contribution is desired, because
an arbitrarily close point may be considered class-defining
under the assumption that the data does not contain erroneous training labels.

Voronoi Boundary Classification: A High-Dimensional Geometric Approach via Weighted Monte Carlo Integration−2 2

If instead we consider ŵ(z) = z −p e−0.5σ z for p ≥ d,this problem can be avoided and the integrals also convergesfor σ = 0. Furthermore, for p = d, the factors of kfj (x)−tk
cancel in the Monte Carlo integration formula of Lemma 2.
Some examples of the classification boundary, includingedge cases, produced by this weight function with different
σ parameters are displayed in the top row of Figure 3. Whenσ is close to 0, the result becomes reminiscent of 1-Nearest
Neighbor classification. This behavior is further discussed
in Section 2. On the contrary, when σ is considerably large,
the role of the distances decreases relative to the role ofthe boundary volumes. Another example that shows a comparison of the decision boundaries between SVM and ouralgorithm for a particular family of data is presented in thebottom row of Figure 3 where we illustrate that our classifier remains largely unchanged even when we introduce an
imbalanced data split between classes. This is in contrast to
the performance of support vector machines and artificial
neural networks (Wang et al., 2016) on imbalanced data.
Optimized implementation We implemented our algorithm with a few optimization considerations in mind.
Firstly, we utilized weight functions with a kfj (m) − tkd
factor to benefit from the previously discussed simplification and to avoid precision problems with computing large
powers of real values in high dimensions. Another set of
optimizations is based on the fact that most of the elementsof the Monte Carlo integration formula can be unfolded andpre-calculated to reduce the computational complexity ofthese elements to constant lookup O(1) during the main
execution of the algorithm. Algorithm 1 presents the complete algorithm, optimized for a test dataset, rather than asingle test point t. This batch implementation allows for
efficient matrix operations as compared to a sequential implementation and was implemented using C++ with the useof the GPU library OpenCL for the most computationallyintensive sections.
In the algorithm, initially we compute all pairwise distances
between test and training data points. Then, during each
iteration of ray sampling, we also precompute all innerproducts between data points and the sampled vector m.
After that, by looking at each pair of points, we can find∗lm(·) in constant time and perform an update to the current
estimate of a corresponding area. By considering the nested
for-loop structure of Alg 1, we observe that the running timeis of order O(N M D + T RD + T N D + T M D + T M N ),where N denotes the number of training data points, M
denotes the number of test data points, D the dimensionof the space containing the data and RD the complexity ofgenerating a sample on the unit sphere in RD .
Approximate Nearest Neighbor Search for Large Data
Note that, while our algorithm consists mainly of vectorized

Algorithm 1 Classification Algorithm
1: Input:
2: Train, test data X ∈ RN ×D , X̂ ∈ RM ×D
3: Train labels Y ∈ {1, . . . , k}N
4: Number of iterations T ∈ N
5: Weight function w : R+ → R+
6: Output:
7: Predicted labels Ŷ ∈ {1, . . . , k}M
8: Compute matrix dst ∈ RN ×M as
9:dstij ← kxi − x̂j k2
10: Initialize matrix areas ∈ RM,k with zeros
11: for it = 1 to T do
12:
Sample random m ∈ RD from U ni(S D−1 )13:Compute vector p ∈ RN as
14:pi ← hm, xi i
15:
Compute vector p̂ ∈ RM as
16:p̂i = hm, x̂i i
17:
for j = 1 to M do // t ← x̂j
18:
Initialize im = 0, lm = +∞19:
for i = 1 to N do
∗
(xi )
20:
lcur ← dstij /(2 · (pi − p̂j )) // lcur ← lm
21:
if 0 < lcur < lm then22:
im ← i, lm ← lcur23:
end if24:
end for25:
if im 6= 0 then26:q ← (pim − p̂j )/sqrt(dstim j ) // q ← hm, ni
27:
Add w(lm ) · lm D−1 /q to areasj,yim
28:
end if
29:
end for
30: end for
31: for j = 1 to M do
32:ŷj ← arg maxl areasj,l33: end forsequential operations and benefits from parallelizability, a
key bottleneck for large datasets is the computation of all
pairwise distances between train and test data points requiring O(M N ) storage and compute time. As an alternative,
we propose an approximate algorithm, which can make useof a chosen approximate nearest-neighbor search algorithm
such as (Muja & Lowe, 2014).
The key difference of the alternative method lies in the wayof finding the face(m), which was introduced in Lemma 2.
For a given sampled direction m we consider the intersection of the ray emanating from a test point t in direction mand use a binary search on the distance to the intersection
from t, with the initial bounds set to 0 and M AX, where thelatter is some chosen maximum value, such as the diameterof the dataset or a value where the weight function’s value
is sufficiently close to zero to be considered insignificant forthe classification task.
Note that, if the intersection occurs at a distance ρ, thenVoronoi Boundary Classification: A High-Dimensional Geometric Approach via Weighted Monte Carlo Integration
2

−2

0.1
w(z) = exp(−0.5z
z2

2 −2

)

1
w(z) = exp(−0.5z
z2

2

)

4

4

4

2

2

2

0

0

0

−2

−2

−2

−4

−4

−2

0

2

−4

4

(a) σ = 0.1, comparable to 1NN.

−4

−2

0

2

−4

4

(b) σ = 1.0.

−4

−2

40

40

40

20

20

20

20

0

0

0

0

−20

−20

−20

−20

−40 −2002040(d) SVM(γ = 0.1) (4 vs 4).

−40

−40
−40 −2002040(e) SVM(γ = 0.1) (8 vs 4).

−40 −20

0

20

40

0

)

2

4

(c) σ = 10.0.

40

−40

−2

10
w(z) = exp(−0.5zz2−40−40 −2002040(f) VBC(σ = 102 ) (8 ”outer” points vs 4 ”inner” points).

Figure 3. Top row: An example of a classification boundary on 6 training points. Bottom row: Comparison of SVM decision boundaries
to our algorithm. The intensity of the color represents the ratio of the largest Voronoi rank divided by the sum of all Voronoi ranks at agiven point (split into ten bins in 3a, 3b, 3c, two bins in 3f (left) and continuous in 3f (right)). Observe that as we double the outer training
data density with a slight displacement from 4 to 8 between figure (d) and (e), the SVM classifier boundary changes drastically. In figure(f) we display the corresponding decision boundary for 8 points for VBC and the full level-set of the decision boundary for 4 points. The
decision boundary in our approach remains largely unchanged since the Voronoi cell geometry is not affected significantly.

points on the ray with distances from the interval [0, ρ] belong to test point’s cell, and points with distances from(ρ, M AX) belong to other cells. Therefore, if t is ourtest point and nearest(·) is the nearest neighbor functionthat returns the closest datapoint to a given one, thenthe binary search is performed over the function h(l) =1[nearest(t + l · m) = t], searching for the argument at
which the function’s value changes from 1 to 0. After theseach is complete, we can decide what cell is getting crossedby the ray by running the neareset neighbor algorithm forthe intersection point one more time.
The nearest neighbor algorithm in this case may be replacedwith an approximate search version, which would allow us
to replace the mentioned part from the algorithm’s complexity with an asymptotically smaller function. The complexity
now becomesO(prep(N, D)+T M log

diam(D)(nn(N, D)+D)) (3)ε0where prep(N, D) is the time needed to construct a data
structure for a nearest-neighbor search, nn(N, D) is thecomplexity of the search over N D-dimensional points(which might be considered as D log N for some algorithms(Arya et al., 1998)) and ε0 is the desired precision for thebinary search.
Convergence to 1-Nearest-Neighbor classifier As wasobserved empirically in Fig. 3, certain weight functions
appear to return classification results comparable to a 1nearest neighbor classifier. Here, we present a theorem
making this observation precise:
Theorem 1. Consider a labeled dataset D={(x1 , c1 ), . . . , (xn , cn )} of data points xi ∈ Rd andcorresponding class labels ci ∈ {1, 2, . . . , k} and atest point t ∈ Rd , where t 6= xi for all i ∈ {1, . . . , n}.
Assume that t only has a single nearest neighbor among{x1 , . . . , xn }. Consider a sequence of weight functions{wn }∞
n=1 , wn : R+ → R+ which are each monotonicallydecreasing and whereZ +∞wn (z) d−2limzdz = 0 for all
0 < z1 < z2n→∞ z
wn (z1 )2(4)Then, for sufficiently large n, the class V BC(t|D, wn ) assigned by the Voronoi Boundary Classifier of t is equalVoronoi Boundary Classification: A High-Dimensional Geometric Approach via Weighted Monte Carlo Integration

to the class 1N N (t|D) assigned by the nearest neighborclassifier for t.
Proof. Please consult the supplementary material.
−0.5z 2 σ −2

Our standard weight function ŵ(z) = e zp
in particular can be considered in the light of the above result. Consider two sequences: wn (z) = ŵ(z|σ = n−1 , p = p0 ≥ 1)
and wn (z) = ŵ(z|σ = σ0 , p = n) for n ∈ N. It is easy
to see that both sequences follow the theorem and hence
converge to nearest neighbor classification when used asweight functions for Voronoi Boundary Classification.

3. Experimental Evaluation
The following weight function was used in all experiments:−22wσ (z) = kzk−d e−0.5σ kzk , where d is the dimensionalityof the data. The parameter σ varies for different datasets.
We noted that in higher-dimensional spaces the classification
performance appeared to improve with a larger σ.
Monte Carlo Convergence Firstly, we investigate theconvergence of our proposed Monte Carlo Integration procedure. For these experiments, we used generated: 2D, 20D
and 1000D normally distributed pointclouds of 100 points
each and also considered the MNIST dataset which is 784
dimensional. We assigned all data points to a single classlabel 1 and study the convergence of the Voronoi Boundary
Rank r1 for a selection of random test points. All tests in thissubsection use σ = 103 as the weight parameter. Each subplot on Figure 4 demonstrates 25 curves of values of r1 (t)as we increase the number of Monte Carlo samples. We
observe after only 105 iterations the approximated integrals
start to stabilize. It should be noted that in this experiment
the points are considered to belong all to 1 class, while ina classification task those calculated areas for each cell are
split among several classes neighboring the test Voronoi
cell. Note that, for MNIST, approximately 104 iterations
were enough for areas to stabilize, while for the synthetic
datasets, 105 iterations were necessary for 20D and 1000D
cases.
Figure 5 furthermore focuses in on the Voronoi Boundary
Rank convergence for two of these datasets, with low- andhigh-dimensional data, but for one randomly picked test data
point and shows the standard deviations for 100 repeatedrandom initializations of the algorithm. As expected fromthe theoretical Monte Carlo integration foundations, thevariance curve can be approximated with a function cn−0.5 ,where n is the number of iterations and this fitting is also
displayed in the figure.
Classification Performance We now study the performance of our algorithm on the following datasets:

1. Anuran Calls – Frog sounds (MFCCs), 22 attributes,
10 species. 7195 data points split into 3598 training,3597 test points (Dheeru & Karra Taniskidou, 2017).
2. MNIST – a classic dataset with 28 × 28 × 1 = 784
dimensional dataset of images with handwritten digits. 60000 training data points, 10000 test data points(Deng, 2012).
3. CIFAR-10 – a dataset with 32 × 32 × 3 = 3072 dimensional image data which corresponds to 10 different
classes; 50000 training data points, 10000 test data
points (Krizhevsky et al., 2014).
Figure 6 first reports results of the convergence of the classification accuracy depending on the number of Monte Carlosamples. As can be seen, after 104 iterations the changes tothe accuracy value become negligible and thus all classification tests were conducted with 104 Monte Carlo samples.
Our algorithm (VBC with 10000 samples) was then tested
against 1NN, (KNN with K ∈ {2, . . . , 10} results are notpresented since examining only one neighbor was always
the most optimal for these datasets), SVM (Chang & Lin,2011) with an RBF kernel and default parameters fromscikit-learn 0.20.0 and RandomForest (Breiman,
2001) with 1500 estimators. The comparison is presented
in Table 1, where results are averaged over 25 runs.
Table 1. Test results

VBC(10k)
1NN
SVM
RF(1500)

F ROGSMNISTCIFAR-10.986±.0002
.982
.903
.974±.0004.969±.0003
.969
.944
.972±.0004.494±.0011
.354
.440
.495±.0016

We observe that Voronoi Boundary Classification obtains
satisfactory results, comparable to the best performance ofthe classical algorithms. In the case of the CIFAR dataset,
where the nearest neighbor algorithm obtains low accuracy,
our algorithm, despite its superficial resemblance to 1NN,
has a significantly better accuracy than 1NN and is in linewith the performance of random forests.
The running time of our algorithm was measured on the
MNIST dataset, on an Intel Core i7-7700HQ processor and
GeForce GTX 1060 graphics card. Performing the classification with 10000 Monte Carlo samples took 4 minutes
53 seconds for all 10000 test images, of which 15 seconds
were spent on initialization of the algorithm and about 0.277
seconds for each ray casting iteration on average.
While these initial tests show the feasibility of VBCclassification on medium-sized datasets, we defer an in
depth performance evaluation to future work since optimizations such as the mentioned approximate nearest neighborsearch implementation and large scale implementation will
require further work.

Voronoi Boundary Classification: A High-Dimensional Geometric Approach via Weighted Monte Carlo Integration
1.6
20

4

0.6

1.4

18
3

0.4

16

1.2
2
0

0.5

14

1

0

0.5

·105

(a) D = 2

1

0

·105

(b) D = 20

0.51(c) D = 1000

0

0.5

·105

1
·104

(d) MNIST

Figure 4. Area convergence results for normally distributed data and MNIST with σ = 103 .
·10−2√
0.2090/ x

2

1.2

std

10.8

0.25

10.7

√
2.4718/ x

0.20

std

10.6

0.15

1.19

1.5

1.18

1

10.5

0.10

1.17

0.5

10.4

0.05

0

10.3

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

×105

0.8

0.00

1

0

0.2

0.4

0.6

0.8

·105

(a) Gaussian; D = 2

1

0

0.2

0.4

0.6

0.8

·105

1
·105

(b) MNIST

Figure 5. Area convergence for a single cell with weight function parameter σ = 103 . The left plot in each subfigure represents theaverage Voronoi Boundary Rank, its standard deviation over 100 executions in grey. The dark gray curve represents a single run andthe blue curve depicts the mean Voronoi Boundary Rank. The second and fourth figure each show the standard deviation and its MSE√
approximation with c/ x.
0.970

0.986

0.500
0.475
0.450
0.425
0.400
0.375
0.350
0.325

0.969

0.984

0.968

0.982

0.967

0.980

0.966
0.965

0.978

0.964
0

0.2

0.4

0.6

0.8

(a) Anural calls

1
·104

0

0.2

0.4

0.6

(b) MNIST

0.8

1
·104

0

0.2

0.4

0.6

0.8

(c) CIFAR-10

1·104

Figure 6. Accuracy plots on the given datasets. We display mean accuracy (blue) and standard deviation (light gray) error bounds over
100 iterations of our algorithm and an example run (dark gray) of our approach.

Performance dependence on training data size In Fig.
7, we study how the accuracy of the algorithms changes withthe size of the training dataset. As all algorithms quickly
improve in accuracy as the training dataset size increases, we
display the difference in performance to the mean accuracy
over all 4 studied algorithms. The results were obtainedby sampling 25 random subsets of the data. The displayed
graphs represent mean performance over those 25 subsets
for each training data size setting. Note that we do notdisplay the results for SVM which were significantly worse
than the displayed values.
The intuition is that our algorithm should benefit, comparedto say 1NN, when it becomes possible to approximate theVoronoi cell geometry.
We note that the Voronoi boundary classification algorithm
performs well on both small and large datasets. On the contrary, most of the other algorithms do not have this featureor require additional tuning, that is dependent on the size ofthe data (as in the case with neural networks). We note that,particularly for CIFAR, and VBC performance for larger σ
parameter, VBC performance is similar for Random Forests
which performed well for large data size, while for small σ
VBC performed similar to 1NN for large training datasetsin line with Theorem 1.

Voronoi Boundary Classification: A High-Dimensional Geometric Approach via Weighted Monte Carlo Integration

knn(n neighbors=1)
vc(sigma=10)·10−2

rf(n estimators=1500)
vc(sigma=100)vc(sigma=1)
vc(sigma=1000)

·10−2

7.5

4

·10−2

5

0.5
2

2.5

0

0

0
−0.5

−2.5
−5

−2

−1

−7.5

−4

−1.5
10

30

100

300

1000

3000

−10
10

30

100

(a) Anural calls

300

1000 3000 10000 30000

10

30

100

(b) MNIST

300

1000 3000 10000 30000(c) CIFAR-10

Figure 7. For each training data set size, we compute the relative difference between the algorithms’ accuracy and mean accuracy over all
4 studied algorithms for that training data set size. Results are furthermore averaged over 25 random training datasets for each training
data size setting.

MNIST
1

2500
number of images

accuracy

0.9

2000

0.8

1500

0.7

1000

0.6

500
0
0

0.2

0.4 0.6
ratio

0.8

1

Figure 8. Correlation between accuracy and confidence value.

Confidence measurement We now consider a measureof confidence in our classifier provided by the formula (S1 −S2 )/(S1 + S2 ), where S1 is the largest Voronoi Boundary
Rank for a given point, and S2 is the second largest Voronoi
Boundary Rank.
The value of this ratio lies in [0, 1], with 1 corresponding
to the fact that only one class was encountered during integration, while 0 indicates that the algorithm has to pickrandomly between at least two classes, due to corresponding integrals being equal. Figure 8 shows how the accuracyof our algorithm correlates with the described confidence
levels on the MNIST dataset. Note that we observe that thealgorithm has accuracy close to 100% for confidence values
above 0.8, as well as a low accuracy of approximately 50%
for confidence values close to 0.
We believe that confidence measures such as the proposed

one may be useful in scenarios when the task of classification is extended by having an option of not assigning a classlabel in case of high uncertainty.

4. Conclusions and Future Work
We have proposed a geometric classification algorithm basedon a weighted integral over Voronoi cell boundaries anddemonstrated the performance of the algorithm on a set ofexample datasts.
We believe that the simplicity of the geometric formulation of this algorithm will allow for extensive analysis interms of error-bounds, behavior depending on training data
distribution assumptions etc. In future work, we hope toevaluate the performance of our approach on larger datasets
and would like to generalize the ray-casting Monte Carlointegration based approach of the proposed algorithm to alarger family of geometric classification algorithms. Finally,
we intend to test our approach also as a hybrid method in
conjunction with Bayesian and Deep Neural Network based
approaches.

5. Supplementary Material and Source CodeThe authors will maintain updated versions of supplementary information for this paper at their respective websites https://people.kth.se/˜vpol and https://people.kth.se/˜fpokorny. Source code will beavailable at https://github.com/vlpolyansky/
voronoi-boundary-classifier

6. Acknowledgements
This work has been supported by the Knut and Alice Wallenberg Foundation.

Voronoi Boundary Classification: A High-Dimensional Geometric Approach via Weighted Monte Carlo Integration

References
Arya, S., Mount, D. M., Netanyahu, N. S., Silverman, R.,
and Wu, A. Y. An optimal algorithm for approximate
nearest neighbor searching fixed dimensions. Journal of
the ACM (JACM), 45(6):891–923, 1998.
Aurenhammer, F. Voronoi diagramsa survey of a fundamental geometric data structure. ACM Computing Surveys(CSUR), 23(3):345–405, 1991.
Biau, G. and Devroye, L. Lectures on the nearest neighbormethod. Springer, 2015.
Breiman, L. Random forests. Machine learning, 45(1):
5–32, 2001.
Chang, C.-C. and Lin, C.-J. Libsvm: a library for support
vector machines. ACM transactions on intelligent systems
and technology (TIST), 2(3):27, 2011.
Davis, P. J. and Rabinowitz, P. Methods of numerical integration. Courier Corporation, 2007.
Deng, L. The mnist database of handwritten digit images
for machine learning research [best of the web]. IEEE
Signal Processing Magazine, 29(6):141–142, 2012.
Dheeru, D. and Karra Taniskidou, E. UCI machine learning
repository, 2017. URL http://archive.ics.uci.
edu/ml.
Döring, M., Györfi, L., and Walk, H. Rate of convergenceof k-nearest-neighbor classification rule. The Journal of
Machine Learning Research, 18(1):8485–8500, 2017.
Doucet, A., De Freitas, N., and Gordon, N. An introduction
to sequential monte carlo methods. In Sequential Monte
Carlo methods in practice, pp. 3–14. Springer, 2001.
Krizhevsky, A., Sutskever, I., and Hinton, G. E. Imagenet
classification with deep convolutional neural networks.
In Advances in neural information processing systems,
pp. 1097–1105, 2012.
Krizhevsky, A., Nair, V., and Hinton, G. The cifar-10 dataset.
online: http://www. cs. toronto. edu/kriz/cifar. html, 2014.
Kurakin, A., Goodfellow, I., and Bengio, S. Adversarial examples in the physical world. arXiv preprint
arXiv:1607.02533, 2016.
Kushilevitz, E., Ostrovsky, R., and Rabani, Y. Efficient
search for approximate nearest neighbor in high dimensional spaces. SIAM Journal on Computing, 30(2):457–
474, 2000.
Muja, M. and Lowe, D. G. Scalable nearest neighbor algorithms for high dimensional data. Pattern Analysis andMachine Intelligence, IEEE Transactions on, 36, 2014.

Sharif Razavian, A., Azizpour, H., Sullivan, J., and Carlsson,
S. Cnn features off-the-shelf: An astounding baseline
for recognition. In The IEEE Conference on Computer
Vision and Pattern Recognition (CVPR) Workshops, June2014.
Sibson, R. A brief description of natural neighbor interpolation (Chapter 2), Interpolating Multivariate Data. John
Wiley, 1981.
Wang, S., Liu, W., Wu, J., Cao, L., Meng, Q., and Kennedy,
P. J. Training deep neural networks on imbalanced data
sets. In Neural Networks (IJCNN), 2016 International
Joint Conference on, pp. 4368–4374. IEEE, 2016.
Weinberger, K. Q. and Saul, L. K. Distance metric learningfor large margin nearest neighbor classification. Journalof Machine Learning Research, 10(Feb):207–244, 2009.