Computing Kernel Density Estimates on the GPU

A histogram counts points into cells and produces a blocky field; a kernel density estimate spreads each point over a neighbourhood with a falloff function and produces a smooth one. The difference matters for the map because a KDE reads as a continuous surface rather than as a grid artefact, and it matters for the GPU because spreading a point over a neighbourhood turns one atomic increment into dozens. This page covers the two ways to compute it — direct scatter with falloff, and a histogram followed by a separable blur — the bandwidth question that decides which is right, and the numerical care the accumulation needs. It is one stage of spatial aggregation in GPU memory.

Two routes to the same density surface Two routes are shown as a four-stage flow. The direct route scatters each point into every cell within the kernel radius, weighting by the falloff, which for a radius of four cells means eighty-one atomic additions per point. The separable route instead bins each point into one cell with a single atomic addition, then blurs the resulting grid horizontally and vertically with a one-dimensional Gaussian, which costs eighteen taps per cell regardless of how many points there are. DENSITY SURFACE · DIRECT vs SEPARABLE Points one record each Bin one atomic per point a blocky grid Blur H 9 taps per cell Blur V 9 more a smooth surface The direct scatter is only competitive when the point count is far below the cell count.
The separable route decouples the cost from the point count: the blur costs the same for four million points as for four thousand, because it operates on the grid rather than on the data.

Runnable reference implementation

The separable-blur route is almost always the right one, and it is two cheap passes rather than one expensive one. First bin points into a grid exactly as in parallel histogram binning; then blur that grid with a one-dimensional Gaussian, horizontally and then vertically.

wgsl
// Pass 2a: horizontal blur. The grid is r32float; taps are precomputed on
// the host so the kernel does no transcendental work per texel.
const TAPS : u32 = 9u;                       // radius 4, odd so it is centred

struct Blur {
  weights : array<vec4<f32>, 3>,             // 9 weights padded to 12 lanes
  dims    : vec2<u32>,
  radius  : u32,
};
@group(0) @binding(0) var<uniform> blur : Blur;
@group(0) @binding(1) var<storage, read>       src : array<f32>;
@group(0) @binding(2) var<storage, read_write> dst : array<f32>;

fn weight_at(i : u32) -> f32 {
  return blur.weights[i / 4u][i % 4u];
}

@compute @workgroup_size(16, 16)
fn blur_h(@builtin(global_invocation_id) gid : vec3<u32>) {
  if (gid.x >= blur.dims.x || gid.y >= blur.dims.y) { return; }
  var total : f32 = 0.0;
  for (var t : u32 = 0u; t < TAPS; t = t + 1u) {
    // Clamp at the edges so the density does not fall off at the grid border.
    let offset = i32(t) - i32(blur.radius);
    let sx = u32(clamp(i32(gid.x) + offset, 0, i32(blur.dims.x) - 1));
    total = total + src[gid.y * blur.dims.x + sx] * weight_at(t);
  }
  dst[gid.y * blur.dims.x + gid.x] = total;
}

The vertical pass is the same kernel with the axes swapped. Separability is what makes this affordable: a 9×9 two-dimensional kernel is 81 taps per cell, while two 9-tap passes are 18 — and the gap widens quadratically with the radius.

typescript
/** Gaussian taps, normalised so the blur preserves the total. */
function gaussianTaps(radius: number, sigma: number): Float32Array {
  const n = radius * 2 + 1;
  const w = new Float32Array(n);
  let sum = 0;
  for (let i = 0; i < n; i++) {
    const x = i - radius;
    w[i] = Math.exp(-(x * x) / (2 * sigma * sigma));
    sum += w[i];
  }
  for (let i = 0; i < n; i++) w[i] /= sum;   // normalise: total density preserved
  return w;
}

Parameter reference

Value Setting here Guidance
Grid format r32float r8unorm saturates at 255 and destroys the tail of the distribution, which is where the interesting structure usually is.
Kernel radius 4 cells A Gaussian is effectively zero past about 3σ; radius ceil(3σ) captures the field without wasted taps.
σ (bandwidth) 1.5 cells Expressed in cells, so it scales with the grid rather than with the map. See the section below.
Normalisation taps sum to 1 Otherwise the blur changes the total, and the legend’s numbers stop meaning anything.
Edge handling clamp Wrapping bleeds density across the antimeridian; zeroing makes the border artificially cold.
Taps per cell, two-dimensional against separable A bar chart of the number of texture or buffer reads per output cell for four kernel radii, comparing a two-dimensional kernel against two separable one-dimensional passes. At radius two the two-dimensional kernel needs 25 taps and the separable pair needs 10. At radius four it is 81 against 18. At radius eight it is 289 against 34. At radius sixteen it is 1089 against 66. The two-dimensional cost grows quadratically while the separable cost grows linearly. TAPS PER CELL Radius 2 25 vs 10 Radius 4 81 vs 18 Radius 8 289 vs 34 Radius 16 1089 vs 66 0 300 600 900 1200 taps 2-D kernel two 1-D passes A Gaussian is separable exactly; most other useful falloff kernels are not, which is why it wins.
Separability is not a small optimisation. At a radius of sixteen the two-dimensional kernel does sixteen times the work for an identical result, which is the difference between a real-time surface and a loading spinner.

Bandwidth in pixels or in metres

The single most consequential choice is what the bandwidth means, and the two answers produce visibly different maps.

Bandwidth in screen pixels keeps the blur the same visual size at every zoom. The density surface looks consistent as the user zooms, which is what most interactive heatmaps do, and it is achieved by keeping σ fixed in grid cells and rebuilding the grid per viewport. The consequence is that the surface is not a property of the data — zooming in reveals structure that was smoothed away, and two screenshots at different zooms are not comparable.

Bandwidth in ground metres makes the surface a property of the data. σ is converted to cells using the current metres-per-cell, so the blur shrinks on screen as the user zooms in and the same features stay visible at the same ground scale. It is what a statistical map needs, because the surface then means something independent of the viewer, and it is what a legend can honestly put a number on.

The rule worth applying is that anything with a legend attached should use ground metres, and anything purely illustrative can use pixels. Mixing them — a pixel bandwidth with a metres-per-cell legend — produces a map whose stated units change as it is zoomed, which is the sort of error that survives review because every individual frame looks reasonable.

Choosing the falloff kernel

The Gaussian is the default for a reason worth understanding, because the alternatives are cheaper and occasionally better.

A Gaussian is the only common kernel that is exactly separable, which is what turns a quadratic cost into a linear one. It is also smooth to every derivative, so the resulting surface has no visible discontinuities at any zoom — no faint rings, no plateaus. Those two properties together make it the right default for anything interactive.

An Epanechnikov kernel — a simple inverted parabola, zero outside its radius — is what statisticians reach for, because it minimises the mean integrated squared error for a given bandwidth. It is not separable, so on a GPU it costs the full two-dimensional tap count, and the theoretical improvement over a Gaussian is small enough that it rarely survives the cost. Where it earns its place is in a server-side computation whose numbers will be published, and where the choice of kernel is something a reader might reasonably question.

A simple box or triangular kernel is cheap and produces visible artefacts: a box kernel gives every point a hard-edged square of influence, which reads as blocky at low densities, and a triangular one leaves faint diamond shapes. Both are separable, so if the surface is purely illustrative and the radius is large, they are defensible — but the visual saving from a Gaussian is usually larger than the arithmetic saving from anything else.

The practical advice is to use a Gaussian on the client and to record which kernel produced any surface whose values are published, because a density value without a stated kernel and bandwidth is not a number a reader can do anything with.

Failure modes

  • The surface saturates to flat white in dense areas. An r8unorm grid clipped at 255. Use r32float and map to colour at the end.
  • The total changes when the blur radius changes. The taps were not normalised. Every legend value is then a function of a tuning parameter, which is worse than being wrong.
  • A cold border around the grid. Edge samples were treated as zero rather than clamped, so cells near the border average in emptiness that is not there.
  • Density leaks across the antimeridian. Edge handling wrapped instead of clamping. In a global grid the two edges are genuinely adjacent, which makes this the one case where wrapping is correct — for a tile grid it never is.
  • The blur is correct but slow. A two-dimensional kernel instead of two one-dimensional passes. At radius 4 that is 81 taps against 18.
What the bandwidth is measured in A comparison of two bandwidth conventions across four properties. A bandwidth fixed in screen pixels keeps the blur visually constant across zoom levels, makes the surface a property of the view rather than of the data, means two screenshots at different zooms are not comparable, and cannot carry an honest legend. A bandwidth fixed in ground metres makes the blur shrink on screen as the user zooms in, makes the surface a property of the data, keeps screenshots comparable, and supports a legend with real units. BANDWIDTH CONVENTION Screen pixels Ground metres Look across zoom constant shrinks in Surface is a property of the view the data Screenshots comparable no yes Legend can state units no yes Anything with a number beside it should use ground metres — no exceptions worth making.
The pixel convention is not wrong, it is illustrative. It becomes wrong the moment a legend is attached, because the units on the legend then change as the reader zooms.

Backend / Python interop note

Where the surface has to be reproducible — published statistics, a report, anything a reader might recompute — the estimate belongs on the server, where scipy provides a reference implementation and the parameters can be recorded alongside the output.

python
import numpy as np
from scipy.ndimage import gaussian_filter

def kde_grid(xs, ys, extent, cells: int, sigma_metres: float) -> np.ndarray:
    """Bin then blur — the same two steps the GPU performs."""
    (x0, y0, x1, y1) = extent
    grid, _, _ = np.histogram2d(xs, ys, bins=cells,
                                range=[[x0, x1], [y0, y1]])
    metres_per_cell = (x1 - x0) / cells
    return gaussian_filter(grid, sigma=sigma_metres / metres_per_cell,
                           mode="nearest")     # "nearest" == clamp at the edges

Two properties are worth carrying in the output’s metadata: the bandwidth in metres, and the cell size. Between them a reader can reconstruct what the surface means, and a client that renders the grid can label its legend without guessing. Recording the kernel alongside them is worth the extra field for the same reason. mode="nearest" matches the clamped edge handling in the shader above, and mismatching the two is a common source of client and server surfaces that agree everywhere except the border.