Reducing the dataset size for KDE

Question

I have GPS data, so 2 coordinates, and I want to estimate the busiest places (i.e. the places with more data points). However, I have a lot of points: currently ~4 million for 12 days, and I will be adding more days with time, and the kernel density estimation (KDE) should be always recomputed when the new data is added.

If I make a random sample of data, would I still get correct KDE? I am a beginner in statistics, but I am afraid that due to Central Limit Theorem it will not be a case.

Or what are other methods of fast KDE computation with the large data sets? Or would I have to discard the historical data and recompute new KDEs only for the recent days, not for all time period?

Answer 1

First off, KDE computation is usually performed with the Fast Fourier Transform (FFT). This algorithm requires $O(N\log N)$ computational time to represent the KDE on a grid of $N$ cells (even in multiple dimensions), making it very fast.

Nevertheless, improvement is possible. The following algorithm allows you to retain the current KDE (which presumably is a raster dataset) and ignore the current data, and then update it with the new data at any time. Thus you could potentially handle far more data than your storage systems permit: you only need to store the KDE and keep updating it, throwing away the data afterwards.

The idea is straightforward: the value of the KDE of a dataset at some point $u$ on your raster is the average of the values that are "spread out" from the data by the kernel function $K.$ We don't even need to consider the formula here; it suffices to know that an observation of the value $z$ at some location $x$ will introduce the value

$$z K(x,u)$$

associated with $u$ and all you have to do is average these values over the data for each $u.$

Fix a $u,$ then, and suppose the KDE of a set of $n$ data values computed for $u$ is $k_0(u)$ and the KDE of a separate set of $m$ data computed for $u$ is $k_1(u).$ Then the KDE for the combined dataset simply is

$$k(u) = \frac{n k_0(u) + m k_1(u)}{n+m}$$

because the numerator reconstructs the sum of all the spread values and the denominator computes the associated average. This operation needs to be done at each point $u$ of the raster.

That leads to the following algorithm. Its input is a raster $K$ of the existing data; the count $n$ of the existing data; the new data and their count $m;$ and any parameters of the kernel density estimate you want to make with the new data. (You can even use a different kernel and bandwidth if you like.)

1. Compute a raster $K^\prime$ representing the KDE of the new data.

2. Output the raster $(nK + mK^\prime)/(n+m).$

The computational cost therefore equals the cost of computing a KDE only for the smaller dataset of new values plus the cost of summing two rasters (which usually is extremely fast). In fact, if the new data values are localized, their KDE may span only a small part of the extent of the original KDE. With the right software you might be able to update just that portion of the original KDE in step (2), thereby saving even more time.