# Fast evaluation of univariate KDE?

I have a KDE obtained with a large sample ($n=10000$). Given that this is defined as the mean of the kernel evaluated at $n$ points, the evaluation of this kernel slows down my code (I need to evaluate it thousands of times). I was wondering if there is a trick to avoid taking the mean of the $n$ points every time I evaluate the KDE?

• can you use an approximate (and random) solution? – user603 Aug 6 '14 at 15:49
• @user603 Yes. Do you mean an approximation in a grid or something like that? – Kako Aug 6 '14 at 15:50
• I have found a naive approximation of my KDE by evaluating a histogram with a large number of bins (500-1000 seem to work well provided my large sample size). So, the only thing I did was to evaluate a step function at every bin. Thanks. – Kako Aug 6 '14 at 16:13
• The naive approximation is formalized and justified by computing the KDE by means of an FFT. To do this, the sample has to be discretized. Since the FFT is so fast, even discretization into millions of bins is not normally an issue. At this point implementation details come to the fore and they suggest that additional information about your situation is essential for choosing a good algorithm. For instance, how do the samples vary from each other? Will they fall into known ranges (allowing bins to be precomputed)? How accurate and precise must the KDE be? Etc, etc. – whuber Aug 6 '14 at 16:31
• If you use R there is a paper by Wickham on KDE that compares different packages: vita.had.co.nz/papers/density-estimation.pdf – Tim Mar 9 '15 at 13:40

In my opinion using a program which evaluates KDE using FFT is quite good choice. This is supported for example in R's density{stats} function or in some dedicated functions in ks package. As was noted before FFT is used a kind of discretization of data is needed (known here as binning, usually implemented as linear binning). But for n like 10000 or so the rounding effect can be totally neglected. KDE is then evaluated almost immediately. You can test for example:

x <- rnorm(10000)
system.time({
# bw can be one of the following: "nrd0", "nrd", "ucv", "bcv", "SJ-ste", "SJ-dpi"
# note that some methods for finding the optimal bw are faster that others
density(x, bw="nrd0")
})

library(ks)
system.time({
# other bandwidth selectors: "hscv", "hlscv", "hpi",  "hns"
kde(x, h=hpi(x))
})