How to speed up Kernel density estimation

Question

Background

A couple of days ago I asked here about: How to reduce number of points for clustering . Instead of reducing number of points, a method "Kernel Density Estimation" (KDE) was suggested to me, it gives right solutions and it is faster then my previous approach.

Question

What I am curious about now is complexity of this method. Maybe I've implemented it in a bad (too naive) way, here is my c++ implementation http://pastebin.com/gtStWjmA (see evalPrivate method). But assuming I have $m$ data points and I want to sample KDE in $n$ points. Then for each evaluation of KDE I have to evaulate kernel function m-times. So my complexity is $O(m^n)$ and that is too much. If my $m=60 000$ and $n = 1000$. Then it tooks ages to sample KDE in order to find its local maximas.

Answer 1

One trick which I used when I implemented KDEs is to limit the effect of kernel to some values.

Suppose you have a sample $x = {x_1, x_2, .. , x_n}$, and some points where you want to estimate the kernel $k = {k_1, k_2, .., k_n}$. Now, without loosing of generality, we can consider $k$ values as being sorted. If not, then you simply sort them and change indexes.

Consider now what a KDE does. Basically, for each $x_i$ you have a probability mass which you want to spread with a symmetric kernel function to the left and to the right. You are interested only to evaluate that probability mass dispersion only onto the values of $k$. Obviously, if you evaluate this kernel function onto all values of $k$, than you will end up with an execution time of $O(|k||x|)$ (not $O(|k|^{|n|})$ as you suggested, consider 20 samples with 20 estimation points is $20^{20}$ = 104 septillion 857 sextillion 600 quintillion operations).

My trick is to not disperse this probability mass to all the points, but only to some of them, the ones which are closer to the given $x_i$ for which the kernel function is evaluated. For some kernel density functions, like uniform, triangular, Epanechnikov this is a natural concept, since for some distance to the left or to the right, the spreaded probability mass equals $0$ so, it does not affect the kernel density estimation at all. For some other kernel functions like Gaussian I established some limits. For example for Gaussian I established the limit as a function of standard deviations. I do not remember the exact value (you can take a look into my code since is publicly available), but outside of a range of some standard deviations to the left of $x_0$ and to the right of it I considered contributed values as being $0$.

Now, the implementation idea is somehow straightforward to follow in practice. I enriched my kernel functions with a minimum value and maximum value (in other words with a distance to the left or to the right of evaluated $x_i$), and using binary search on vector $k$, I find the range of values from $k$ which has to be updated.

Obviously this can produce errors. But, if you have enough points the effect I found is not a problem in practice. And if you implement it well, using range around $x_i$ as a parameter (I did not do it, but I might if I will need that), you will have a parameter which can be use a a compromise between the precision of the evaluation of KDE and speed. Note also that for some kernel functions this parameter has no effect, because the KDE is precise.

Answer 2

Suppose you have $x_1,\dots,x_m$ data points, called source points hereafter, and $z_1,\dots,z_n$ points you want to estimate at, called query points hereafter. A naive implementation of

$$ \hat f(z_j)=\frac 1{mh}\sum_{i=1}^mK\left(\frac{x_i-z_j}{h}\right)$$

is then $\mathcal{O}(mn)$ if you want to evaluate it for all $n$ query points.

Options in 1D

In 1D then sort both the source and query points in $\mathcal{O}(m\log m)$ and $\mathcal{O}(n\log n)$ time. Then pass through both simultaneously and do not evaluate pairs where $K\left((x_i-z_j)/h\right)<\epsilon$ for some small $\epsilon$. For the latter you do not need to check a lot of values as the data is sorted and you pass from left to right. Thus, this step is $\mathcal O(m)$ or $\mathcal O(n)$ depending on which is largest. This is even exact in some cases where the kernel has finite support (like the Epanechnikov kernel).

An even faster option if you have substantial amount of data is to use binning.

Options in > 1D

One option in higher dimensions is the dual-tree method like the one suggested by

Gray, Alexander G., and Andrew W. Moore. "Rapid Evaluation of Multiple Density Models." In AISTATS. 2003.

On example is to use a k-d tree for both the source and query points. Then you travers down both trees simultaneously exploiting that:

• $K\left(\frac{x_i-z_j}{h}\right)$ is virtually zero for points in the current node of source points and query points. Thus, you can use an approximation when it is below a certain threshold.
• the distance between source points are small in some nodes in the source point tree. Thus, $K\left(\frac{x_i-z_j}{h}\right)$ is almost a constant for a node with query points sufficiently far away. Again, you can use an approximation when the distance is sufficiently small in the source node and query node is far away.

The above may not even introduce errors if the kernel has finite support. The algorithm will have $\mathcal{O}(m\log m)$ and $\mathcal{O}(n\log n)$ run-times (depending on which is largest). I have an example here of very similar problem that occurs in particle smoothers. Here is an example of run-times of the naive method versus the dual-tree method

##          method
## N         Dual-tree      Naive Dual-tree 1
##   384       0.01341  0.0008045     0.00529
##   768       0.04257  0.0028020     0.01418
##   1536      0.07590  0.0103289     0.02291
##   3072      0.08270  0.0395034     0.03428
##   6144      0.07206  0.1071989     0.05398
##   12288     0.03143  0.4850711     0.09752
##   24576     0.05360  1.7199726     0.18629
##   49152     0.10605  6.9812358     0.36645
##   98304     0.19286 31.1321333     0.77592
##   196608    0.39916         NA          NA
##   393216    0.86169         NA          NA
##   786432    1.80240         NA          NA
##   1572864   3.66839         NA          NA
##   3145728   8.10332         NA          NA

N in the above is $N=m=n$ and the run-times are in seconds. The last column is a single threaded version.