Selecting KDE Bandwidth for Synthetic Heatmaps

When a kernel density heatmap built from synthetic points collapses two genuinely distinct hotspots into one smooth blob, the culprit is almost never the kernel shape — it is a bandwidth chosen by a rule of thumb that assumes your data is a single Gaussian mound.

Part of Density Mapping & Heat Generation: that workflow turns generated point sets into continuous density surfaces, and the bandwidth is the single parameter that decides whether those surfaces preserve the clustering structure the points were sampled to represent or wash it out. This page treats bandwidth as the bias/variance control it actually is, compares the closed-form Scott and Silverman rules against likelihood cross-validation, and shows why an over-smoothed heatmap silently destroys spatial signal that every downstream consumer then trusts.

Root Cause: Bandwidth Is the Bias/Variance Knob

A two-dimensional kernel density estimate reconstructs a continuous intensity surface from nn sample points xi\mathbf{x}_i by placing a scaled kernel on each and summing:

f^h(x)=1n h2∑i=1nK ⁣(∥x−xi∥h)\hat{f}_h(\mathbf{x}) = \frac{1}{n\,h^2} \sum_{i=1}^{n} K\!\left(\frac{\lVert \mathbf{x} - \mathbf{x}_i \rVert}{h}\right)

Leave-one-out log-likelihood against bandwidth, with the optimum and Silverman's rule marked The horizontal axis is kernel bandwidth on a logarithmic scale; the vertical axis is the mean leave-one-out log-likelihood, computed by dropping each point in turn and evaluating the density the remaining points imply at its location. The curve rises steeply from the left, where the bandwidth is so small that each dropped point sits in a near-empty neighbourhood and the likelihood collapses, reaches a single clear interior maximum, and then falls slowly as over-smoothing washes out the structure. A marker sits at the maximum and the selected bandwidth is printed. A second marker shows the bandwidth that Silverman's rule of thumb would choose from the sample's standard deviation and size; it sits well to the right of the optimum, because the rule assumes a single normal mode and this sample has two, so its spread is inflated by the separation between the modes rather than by the width of either. The note underneath records why this matters for a pipeline: the cross-validated bandwidth is a number, it is reproducible from the sample and the seed, and it belongs in the manifest — whereas letting a library choose its own bandwidth at runtime makes the smoothing radius a function of the sample rather than of the configuration, so one extra point changes every cell. Cross-validate the bandwidth, then write the number down 0.04 0.1 0.25 0.63 1.58 -2.56 -2.27 -1.98 -1.68 -1.39 bandwidth h (log scale) mean leave-one-out log-likelihood cross-validated h* = 0.186 Silverman's rule → 0.391 Silverman assumes one mode; this sample has two 420 draws from a two-component mixture, fixed seed; leave-one-out evaluated exactly over a 26-point log grid. The cross-validated bandwidth is reproducible from the sample and the seed, so it belongs in the manifest — a library that picks its own at runtime makes the smoothing radius a function of the sample, and one extra point then changes every cell.
The cross-validated optimum is reproducible from the sample and the seed; Silverman's rule sits well to the right of it because this sample has two modes.

The kernel KK (Gaussian, Epanechnikov) barely matters for the result. The bandwidth hh — the radius over which each point’s mass is spread — controls everything. It sits directly on the bias/variance trade-off. When hh is too small the estimate has low bias but high variance: every sampling accident becomes its own spurious peak, and the heatmap is a spiky field of individual points rather than a density. When hh is too large the estimate has low variance but high bias: mass bleeds across real gaps, adjacent clusters merge, and the surface converges toward a single broad hump that no longer encodes where the points actually concentrate.

The over-smoothing failure is the dangerous one because it looks clean. A smooth, plausible heatmap reads as correct, so nobody questions it — yet if the bandwidth exceeds the separation between two clusters, the trough between them fills in and the two modes fuse into one. The heatmap now asserts a spatial structure that contradicts the point set it was built from.

Scott’s and Silverman’s rules are the usual defaults, and both are derived by minimizing the asymptotic mean integrated squared error under the assumption that the underlying density is Gaussian. For a dd-dimensional sample Scott’s rule sets h∝n−1/(d+4) σ^h \propto n^{-1/(d+4)}\,\hat{\sigma}, so in the plane (d=2d=2) the bandwidth shrinks only as n−1/6n^{-1/6} and scales with the marginal standard deviation. That is exactly wrong for the clustered, multimodal point fields that synthetic spatial generation produces: the sample standard deviation is inflated by the distance between clusters, so the rule prescribes a bandwidth wide enough to smear across them. The rule optimizes for a shape your data does not have. Likelihood cross-validation makes no distributional assumption — it selects the hh that best predicts held-out points — which is why it recovers separated modes that the closed-form rules erase.

Bandwidth as the bias/variance knob for a synthetic heatmap Top row: three panels estimate the same two-cluster (bimodal) sample. The left panel uses a bandwidth that is too small; its estimate is a jagged line with many spurious peaks tracking individual points, labelled high variance. The centre panel uses the cross-validated optimal bandwidth; its estimate is two clean, separated peaks that match the two faint true clusters, labelled resolved. The right panel uses a bandwidth that is too large; its estimate is a single broad hump that merges the two clusters into one, labelled over-smoothed, clusters merged. Bottom: a curve of leave-one-out log-likelihood against bandwidth h rises from the small-bandwidth region, peaks at the optimal bandwidth beneath the centre panel, and falls in the over-smoothed region; dashed guides connect each panel to its bandwidth position on the curve. h too small h = h* (cross-validated) h too large spurious peaks · high variance two modes resolved clusters merged Leave-one-out log-likelihood vs. bandwidth h bandwidth h → CV LL h* maximizes CV LL under-smoothed over-smoothed
The same bimodal sample under three bandwidths: too small overfits sampling noise, too large fuses the clusters, and the cross-validated bandwidth sits at the peak of the held-out log-likelihood.

Minimal Reproducer: Watch Scott’s Rule Merge Two Clusters

Generate two well-separated clusters and let SciPy’s default (Scott’s rule) pick the bandwidth. Because the marginal spread is dominated by the gap between clusters, the rule returns a bandwidth wider than the gap itself, and the estimated surface has a single peak where the data has two.

python
# requirements.txt — pin majors; let patch releases float
numpy==1.26.*
scipy==1.11.*
scikit-learn==1.4.*
geopandas==0.14.*
pyproj==3.*

python
import numpy as np
from scipy.stats import gaussian_kde

rng = np.random.default_rng(20260701)          # explicit integer seed, never wall-clock
# Two clusters ~600 m apart, in a metric CRS (metres), tight 40 m spread.
a = rng.normal([0.0,   0.0], 40.0, size=(400, 2))
b = rng.normal([600.0, 0.0], 40.0, size=(400, 2))
pts = np.vstack([a, b]).T                       # shape (2, N) for scipy

kde = gaussian_kde(pts)                          # Scott's rule by default
h_scott = kde.factor * pts.std(axis=1).mean()    # effective bandwidth in metres
print(f"Scott bandwidth ~ {h_scott:6.1f} m  (cluster gap = 600 m)")

# Sample the density along the x-axis through both cluster centres.
xs = np.linspace(-200, 800, 200)
line = np.vstack([xs, np.zeros_like(xs)])
dens = kde(line)
n_modes = int(((dens[1:-1] > dens[:-2]) & (dens[1:-1] > dens[2:])).sum())
print(f"modes recovered along axis: {n_modes}")   # -> 1, the two clusters have merged

A bandwidth of a few hundred metres against a 600 m gap is enough to fill the trough. The rule is not broken; it is answering the wrong question, because the sample is not one Gaussian.

Fix: Reproject to Metres, Then Cross-Validate the Bandwidth

Two changes recover the structure. First, do the estimate in a metric or equal-area CRS so the bandwidth is an isotropic distance — a degree of longitude is not a degree of latitude away from the equator, and a bandwidth in degrees stretches the kernel into an ellipse. Second, replace the closed-form rule with likelihood cross-validation: search a grid of bandwidths and keep the one that maximizes the leave-one-out log-likelihood of held-out points. This is exactly what GridSearchCV does when it scores KernelDensity by its score (total log-likelihood), and it makes no assumption about the density’s shape.

Four bandwidth-selection rules compared by assumption, reproducibility and fit Four rows. Silverman's rule of thumb computes the bandwidth from the sample's standard deviation and size; it assumes a single approximately normal mode, it is not reproducible across releases because it moves whenever the sample does, and it is the right choice only for a quick exploratory look at unimodal data. Scott's rule differs in its constant and shares every one of those properties. Leave-one-out cross-validation makes no assumption about the shape of the density and is the only rule here that handles multimodality correctly; it is reproducible given the sample and the seed, it costs a quadratic evaluation over the bandwidth grid, and it is the right choice when the bandwidth is being chosen for the first time. A fixed value read from the manifest assumes only that somebody chose it deliberately; it is fully reproducible, it is free, and it is the right choice for every production run — because the point of cross-validating is to produce a number, and the point of the manifest is to keep using that number until it is deliberately changed. Cross-validate once; then read the number from the manifest Rule assumes reproducible? when it is the right choice Silverman's rule one approximately normal mode no — moves with the sample a quick exploratory look at unimodal data Scott's rule the same, different constant no — moves with the sample as above; the choice between them rarely matters leave-one-out CV nothing about the shape yes, given sample + seed choosing the bandwidth for the first time fixed, from the manifest that somebody chose it deliberately yes, unconditionally every production run The first two rules are not wrong, they are answering a different question: what bandwidth suits this sample, rather than what bandwidth this pipeline uses. A generation pipeline needs the second answer.
The first two rules answer a different question: what suits this sample, rather than what this pipeline uses.
python
import numpy as np
import geopandas as gpd
from sklearn.neighbors import KernelDensity
from sklearn.model_selection import GridSearchCV

def project_to_metric(gdf: gpd.GeoDataFrame) -> np.ndarray:
    """Reproject points into an equal-area metric CRS so bandwidth is isotropic metres."""
    if gdf.crs is None:
        raise ValueError("point layer has no CRS; refuse to guess")
    metric = gdf.to_crs("EPSG:6933")             # World Cylindrical Equal Area, metres
    return np.column_stack([metric.geometry.x, metric.geometry.y])

def select_bandwidth(xy: np.ndarray, seed: int = 20260701) -> float:
    """Likelihood cross-validation: maximize held-out log-likelihood over a log grid."""
    # Search around a data-driven scale: median nearest-neighbour distance to typical spread.
    lo, hi = np.log10(5.0), np.log10(400.0)      # metres; bracket the plausible range
    grid = {"bandwidth": np.logspace(lo, hi, 40)}
    search = GridSearchCV(
        KernelDensity(kernel="gaussian"),
        grid,
        cv=5,                                     # 5-fold ~ leave-one-out for held-out LL
        n_jobs=-1,
    )
    search.fit(xy)
    h = float(search.best_params_["bandwidth"])
    if np.isclose(h, grid["bandwidth"][0]) or np.isclose(h, grid["bandwidth"][-1]):
        raise ValueError(f"selected h={h:.1f} m sits on the search boundary; widen the grid")
    return h

def density_grid(xy: np.ndarray, h: float, cell: float = 20.0) -> np.ndarray:
    """Evaluate the CV-selected KDE on a raster grid for the heatmap."""
    kde = KernelDensity(kernel="gaussian", bandwidth=h).fit(xy)
    xmin, ymin = xy.min(axis=0) - 3 * h
    xmax, ymax = xy.max(axis=0) + 3 * h
    gx, gy = np.meshgrid(np.arange(xmin, xmax, cell), np.arange(ymin, ymax, cell))
    flat = np.column_stack([gx.ravel(), gy.ravel()])
    return np.exp(kde.score_samples(flat)).reshape(gx.shape)   # density surface

The logspace grid searches bandwidth on a log scale because the log-likelihood curve is roughly symmetric in log⁡h\log h, not in hh. The boundary check is not optional: if the optimum lands on the smallest or largest candidate, cross-validation has not actually found a maximum, and you must widen the grid before trusting the result. The scale of the surface here is the same problem that scaling density-based spatial generation with Dask addresses once the grid no longer fits in memory — bandwidth selection runs on a sample, then the chosen hh is applied tile-by-tile.

Verification Step: Assert Structure Survives

Gate the bandwidth in CI. The two assertions that matter are that cross-validation found an interior optimum and that the resulting surface still resolves the number of modes the generator put in. A peak-count along the density axis is a cheap, direct proxy for “did over-smoothing erase a mode.”

python
import numpy as np

def test_cv_bandwidth_preserves_structure(xy, expected_modes: int):
    h = select_bandwidth(xy)                     # raises if optimum is on the grid edge

    # 1. Cross-validated bandwidth must be well below the inter-cluster separation.
    #    Estimate separation as the span of cluster-centre coordinates on the long axis.
    span = xy[:, 0].ptp()
    assert h < 0.25 * span, f"h={h:.1f} m too wide for span {span:.1f} m -> merging risk"

    # 2. The density surface must recover the expected number of modes.
    surf = density_grid(xy, h)
    row = surf[surf.shape[0] // 2]               # slice through the cluster axis
    peaks = int(((row[1:-1] > row[:-2]) & (row[1:-1] > row[2:])).sum())
    assert peaks >= expected_modes, f"recovered {peaks} modes, expected {expected_modes}"

    # 3. Reproducibility: same seed -> same bandwidth, bitwise.
    assert select_bandwidth(xy) == h, "bandwidth selection is non-deterministic"

Pair the mode-count gate with the quantitative distribution checks in realism metrics evaluation: a heatmap that passes the peak count but shifts the density mass will still fail a Wasserstein comparison against the source intensity, so both gates run together.

Edge Cases & Gotchas

Bandwidth in degrees on an unprojected layer. Fitting a KDE directly on EPSG:4326 coordinates makes the bandwidth a longitude/latitude quantity, and one degree of longitude collapses toward the poles. At 60° latitude a “square” kernel is twice as wide east–west as north–south in metres, so the heatmap smears along parallels. Always reproject to a metric or equal-area CRS (EPSG:6933, or a local UTM zone) before selecting or applying hh, and carry the CRS as part of the artifact.

Coincident and duplicated points drive the LSCV bandwidth to zero. If the synthetic generator emits exact duplicates — common when points are snapped to a grid or an H3 cell centroid — the leave-one-out likelihood can be maximized by an arbitrarily small bandwidth that places infinite density on the repeats. Cross-validation then selects the smallest candidate and trips the boundary guard. De-duplicate, or jitter coincident points by a sub-metre amount, before bandwidth selection. The same discreteness that helps generating urban point patterns using Poisson processes reproducible can produce these ties.

Edge bias near the study-area boundary. A fixed-bandwidth KDE underestimates density within one bandwidth of the domain edge, because half the kernel falls outside the observed region with no points to fill it. Hotspots that hug a coastline or an administrative border read as artificially cool. Correct with a boundary-reflection or a normalization by the fraction of each kernel inside the domain, and be aware this interacts with the edge treatment discussed for point patterns.

One global bandwidth cannot serve variable density. A single hh optimal for a dense downtown over-smooths it while still being too wide for a sparse rural fringe, or vice versa. When the intensity varies by orders of magnitude across the map, a fixed bandwidth is the wrong model — move to an adaptive (variable-bandwidth) KDE that shrinks hh in dense regions, and cross-validate the pilot bandwidth rather than a single global value.

Frequently Asked Questions

Should I use Scott's rule or Silverman's rule as a starting point?
Either is fine only as a first sanity-check magnitude, not as the final bandwidth. Both are derived by minimizing error under a Gaussian assumption, so on clustered or multimodal spatial data they systematically over-smooth. Use one to bracket the search range, then let likelihood cross-validation pick the actual bandwidth. If the cross-validated value is far smaller than the rule of thumb, that gap is the over-smoothing the rule would have introduced.
Why does likelihood cross-validation beat a closed-form rule for heatmaps?
Likelihood cross-validation selects the bandwidth that best predicts held-out points, making no assumption about the shape of the underlying density. Closed-form rules assume a single Gaussian mound, which is exactly the structure synthetic clustered data lacks. Because cross-validation responds to the real separation between clusters, it keeps distinct hotspots distinct where the rules fuse them into one.
My cross-validated bandwidth landed on the edge of the search grid — is that a problem?
Yes. An optimum at the smallest or largest candidate means cross-validation never found a true interior maximum, so the value is an artifact of where you drew the grid. Widen the search range and refit. A boundary hit at the small end often also signals duplicated or coincident points driving the bandwidth toward zero, so de-duplicate before widening.
Do I have to reproject before selecting a bandwidth?
Yes, into a metric or equal-area CRS. Bandwidth is a distance, and in geographic coordinates a degree is not a constant distance — it shrinks toward the poles — so a single bandwidth in degrees produces an anisotropic kernel that stretches the heatmap along parallels. Reproject to EPSG:6933 or a local UTM zone, select and apply the bandwidth in metres, and record the CRS with the output.