Kernel Density Estimation and the Bandwidth Trade-off

A histogram approximates a density by counting points in fixed cells, and its shape depends on where the cell edges happen to fall. Kernel density estimation removes the edge placement by placing a small bump at every observation and summing the bumps, which yields a smooth estimate defined everywhere. The cost is one tuning parameter that controls how wide the bumps are, and the way that parameter must shrink as the sample grows is the whole story of the method’s strengths and its failure in high dimensions.

One kernel per observation

For n points in d dimensions and a kernel function K that is non-negative and integrates to one, the estimate at a query point is

$$\hat{f}(x) = \frac{1}{n h^d} \sum_{i=1}^{n} K!\left(\frac{x - x_i}{h}\right)$$

The factor h^d keeps the total mass equal to one as the bandwidth changes; without it the estimate would integrate to h^d. The query point contributes to the sum through every observation, so evaluating the estimate at m points costs n m kernel evaluations unless an index or a grid approximation is used, and the estimate itself is a density: it is non-negative and integrates to one over the domain.

A kernel that is a probability density makes the estimate a mixture distribution with n components, one per observation, all sharing the same shape and bandwidth. That is precisely the object compared with a fitted mixture model in Gaussian Mixtures and the EM Lower Bound: the mixture chooses a small number of components and fits their positions, shapes, and weights, while the kernel estimator fixes one component per observation and gives up all adaptivity except the common bandwidth. In low dimensions that is a good trade, because a non-parametric estimate cannot mis-specify the shape of the distribution.

For more than one dimension, the kernel is usually a product of one-dimensional kernels, K(u) = prod_j k(u_j), applied after scaling each coordinate by its own bandwidth. The product form assumes the density is not rotated: a strongly correlated pair of features produces a diagonal kernel that never matches the ridge in the data. A full bandwidth matrix can represent rotation and differing spreads, but it has d(d+1)/2 parameters, so it needs a sample size well above d^2 before the estimate is stable, and most practical use stays diagonal.

Bias, variance, and the rate at which the bandwidth shrinks

The estimate has a bias and a variance that move in opposite directions with h. For a symmetric kernel with second moment sigma_K^2 = integral u^2 K(u) du and roughness R(K) = integral K(u)^2 du, a Taylor expansion around a point where the density is smooth gives

• bias^2 approximately (h^4 / 4) sigma_K^4 f''(x)^2
• variance approximately f(x) R(K) / (n h)

so the mean squared error is a sum of a term that grows with h and a term that falls as 1/(n h). A large bandwidth averages over a wide neighbourhood: the estimate becomes smooth, low variance, and biased toward a local average, flattening peaks and filling valleys. A small bandwidth tracks the sample: low bias, but each bump stands alone and the estimate is spiky, with visible individual observations.

Minimizing the sum gives the optimal rate h proportional to n^{-1/5} in one dimension. The exponent is the important part. Halving the bandwidth requires 2^5 = 32 times as many observations, and improving the estimate by an order of magnitude in error requires a four-order-of-magnitude increase in sample size, because the MSE falls as n^{-4/5}. Kernel estimators converge slowly compared with parametric estimators, which is the price of not assuming a functional form.

Choosing the constant is the practical difficulty. A rule of thumb for a Gaussian kernel substitutes the normal density’s curvature into the optimum, giving h = 0.9 min(sigma_hat, IQR/1.34) n^{-1/5} in one dimension, where the interquartile range divided by 1.34 estimates the scale robustly so that a few outliers do not inflate the bandwidth. The corresponding multivariate rule uses n^{-1/(d+4)} times a per-coordinate scale. Both assume the underlying density is roughly normal, so they over-smooth multimodal data, whose true optimum is driven by the curvature of the modes rather than by the overall spread.

The kernel function itself matters much less. Among non-negative kernels with unit variance, the Epanechnikov kernel minimizes the asymptotic mean squared error, and the Gaussian kernel loses only a few percent of efficiency relative to it. The bandwidth dominates the error by a wide margin, which is why a substantial literature exists on choosing it and almost none on choosing the shape. Cross-validation on the estimated likelihood is one principled route: maximize the leave-one-out likelihood with respect to h, which estimates the density with each point excluded and then evaluates it at that point, so a bandwidth that oversmooths is penalized by the poor likelihood of the points it smoothed over.

The dimensional penalty

The optimal bandwidth in d dimensions scales as n^{-1/(d+4)}, and the mean squared error as n^{-4/(d+4)}. Both degrade quickly. To halve the bandwidth in ten dimensions requires 2^{14} = 16384 times as many observations; in twenty dimensions the factor is 2^{24}, about seventeen million. Even before that, the ratio is nearly insensitive to the sample size: at n = 1000 and d = 10, the factor n^{-1/14} is approximately 0.61, so an optimal bandwidth is already six tenths of the coordinate scale, and doubling the sample multiplies it by 2^{-1/14}, a change of about five percent. The estimate has effectively degenerated into one broad bump fitted to the overall location of the data, with all the internal structure smoothed away.

The cause is geometric rather than statistical in origin: the volume of a ball grows as r^d, so a fixed number of neighbours is spread over an exponentially larger region as dimensionality rises, and no bandwidth can be simultaneously local and populated. The practical consequences are that a likelihood threshold for detecting unusual points loses its power, because the estimated density varies little across the space, and that feature selection or a projection to fewer dimensions becomes a precondition rather than an optimization. The earlier treatment of distance concentration in What a Distance Function Decides About Cluster Shape is the same phenomenon seen from the distance side.

import numpy as np

def kde_scores(X, queries, factor=None):
    n, d = X.shape
    if factor is None:                      # Scott's normal-reference factor
        factor = n ** (-1.0 / (d + 4.0))
    h = factor * X.std(0)                   # per-coordinate bandwidths
    out = np.empty(len(queries))
    for i, q in enumerate(queries):
        u = (q - X) / h                     # standardize each coordinate
        k = np.exp(-0.5 * (u ** 2).sum(1)) / np.prod(h)   # product Gaussian kernel
        out[i] = k.mean()                   # pointwise density estimate
    return out, np.log(np.maximum(out, 1e-300))

Boundary bias and the shape of the domain

A kernel centred near the edge of the support extends outside it, where there can be no observations. The mass that falls beyond the boundary is simply missing from the estimate, so the density is systematically underestimated in a strip of width comparable to h along the edge. The magnitude is not a small refinement: with an interior bias of order h^2, the bias at a boundary with a second-order kernel is of order h, which dominates for the modest bandwidths used in practice. Increasing the sample does not remove it, because h shrinks slowly; the estimate converges to the wrong value at the edge.

Three repairs are common. Reflection mirrors the observations across the boundary, so that the kernels reaching outside are cancelled by mirrored kernels reaching in, and restores the correct asymptotic order at straight edges. A boundary kernel modifies the shape near the edge so that it stays inside the support, at the cost of no longer being a probability density and of needing a different kernel for each distance from the boundary. Local linear estimation fits a line, or a plane, by weighted least squares at each query point, which removes the leading bias term everywhere and is the standard recommendation for data with curved or irregular boundaries. For strictly positive variables, transforming to a logarithmic scale makes the boundary infinite and removes the problem altogether, at the cost of estimating a density in the transformed coordinates, as the earlier discussion of scaling and metrics in distance and cluster shape notes.

Modes, outliers, and why local methods replace it

Two of the downstream uses of a density estimate follow directly from its arithmetic. The modes, which are the local maxima of the estimate, are the centres used by mean-shift clustering: a gradient step on the estimate moves each point toward the weighted mean of its neighbours, and the basins of attraction of the modes define the clusters. The bandwidth is then the entire model, controlling both smoothing and the number of clusters found, which is the same parameter tension visible in DBSCAN between the radius and the cluster count. Computing the modes requires evaluating the density at each iteration, so the cost is proportional to n^2 per step unless an index or a binned approximation is used.

The second use is scoring: a point whose estimated density is far below the rest is an outlier, and thresholding the log-density gives a ranking. In low dimensions that ranking is well behaved and has a clear probabilistic interpretation. In high dimensions the bandwidth problem above makes it useless, because the density differences between typical and atypical points are compressed. The replacements estimate local density from a fixed number of neighbours instead of a fixed radius, taking the ratio of a point’s average distance to its k nearest neighbours and the corresponding averages for those neighbours, so that a point in a sparse region is compared with its own neighbourhood rather than with the whole dataset, or use randomized splitting rules whose path lengths are short for isolated points. Both retain the adaptability that a single global bandwidth cannot provide, and both introduce their own parameters, so the trade is not one of tuning burden but of which local scale the score is computed over.

Kernel estimation is not the only smooth non-parametric alternative to a histogram; a histogram is the special case of a box kernel placed on a grid, with bin width chosen at the slower rate n^{-1/3} and an appearance that changes with the origin of the grid. The comparison to the embedding methods is instructive for a different reason: where a density estimate is a well-defined object with a convergence rate but a fixed neighbourhood, an embedding is a picture with no rate at all and a neighbourhood chosen per point. Neither statement says which is the right tool, only what each one can be trusted to report.