Mean Shift: Complete Guide — Principles, Examples, and Python Implementation

Summary

Mean Shift is a non-parametric clustering algorithm based on kernel density estimation. Unlike K-Means, it does not require the user to set the number of clusters in advance. Instead, Mean Shift automatically searches for the modes (peaks) of the underlying data density distribution, and assigns each observation point to the mode toward which it converges. This density-based approach allows it to detect clusters of arbitrary shape and adapt the number of clusters to the real structure of the data.

In this complete guide, we explore the mathematical principles of Mean Shift, its geometric intuition, its practical implementation with scikit-learn, hyperparameter tuning, and four concrete use cases.

Mathematical Principle of Mean Shift

Kernel Density Estimation

The theoretical foundation of Mean Shift rests on Kernel Density Estimation (KDE). Given a set of n data points {x₁, x₂, …, xₙ} in ℝᴰ, the kernel density estimator at point x is written:

f̂(x) = (1 / n·hᴰ) · Σᵢ₊₁ⁿ K((x – xᵢ) / h)

where:

• K is the kernel function (usually a Gaussian or an Epanechnikov kernel),
• h is the smoothing parameter called bandwidth,
• d is the dimension of the data space.

The KDE assigns each point x a density value proportional to the number of observation points located in its neighborhood, weighted by the kernel. The more a point lies in a dense data region, the higher its estimated density.

Gradient ascent toward the local mode

The goal of Mean Shift is to find the modes of the estimated density function f̂(x). A mode is a local maximum of the density function — that is, a point where the gradient vanishes and the Hessian matrix is negative definite.

To reach a mode, Mean Shift uses gradient ascent on the density estimator. The gradient of f̂(x) can be expressed as:

∇f̂(x) = (2 / n·hᴰ⁺²) · Σᵢ₊₁ⁿ (x – xᵢ) · g(||(x – xᵢ)/h||²)

where g is the derivative function of the kernel (profile derivative). This gradient always points in the direction of the steepest increase in local density.

The mean shift vector

The mean shift vector mₕ(x) at point x is defined as the weighted average of neighbors within the bandwidth h window, recentered relative to x:

mₕ(x) = [Σᵢ₊₁ⁿ xᵢ · g(||(x – xᵢ)/h||²)] / [Σᵢ₊₁ⁿ g(||(x – xᵢ)/h||²)] – x

This vector has a fundamental property: it is proportional to the normalized gradient of the density estimator:

∇f̂(x) ∝ f̂_G(x) · mₕ(x)

where f̂_G(x) is the KDE using the kernel G(z) = g(||z||²). This means the mean shift vector always points in the direction of increasing density.

Convergence toward density modes

The Mean Shift algorithm proceeds by successive iterations. At each step t, the position of a point is updated according to:

x⁽ᵗ⁺¹⁾ = x⁽ᵗ⁾ + mₕ(x⁽ᵗ⁾)

that is:

x⁽ᵗ⁺¹⁾ = [Σᵢ₊₁ⁿ xᵢ · g(||(x⁽ᵗ⁾ – xᵢ)/h||²)] / [Σᵢ₊₁ⁿ g(||(x⁽ᵗ⁾ – xᵢ)/h||²)]

This iteration is guaranteed to converge to a local mode of the estimated density, under general conditions on the kernel. The convergence proof relies on the fact that the sequence {f̂(x⁽ᵗ⁾)} is monotonically increasing and bounded above.

Once all points have converged to their respective modes, all points that converge to the same mode are grouped into a single cluster. The number of clusters is therefore determined automatically: it is the number of distinct modes found.

Geometric intuition: marbles in a valley

To understand Mean Shift without formulas, imagine the following scenario.

You have a three-dimensional terrain whose altitude represents the density of your data. Regions where points are numerous form hills; empty regions are valleys.

Now, place a marble at each data point, as if you were scattering marbles across this landscape. Each marble rolls upward along the steepest slope (this is the equivalent of gradient ascent). Neighboring marbles naturally converge toward the same hilltop.

All marbles that end up at the same summit belong to the same cluster. Points on the sides of the hill slide toward the local center of gravity, then continue upward until they reach the peak.

This analogy explains several key properties:

• No need to specify the number of clusters: the number of natural summits automatically determines the number of groups.
• Clusters of arbitrary shape: unlike K-Means which favors spherical clusters, Mean Shift’s basins of attraction can follow complex shapes.
• Bandwidth sensitivity: if the band is too wide, all marbles converge to a single summit. If it is too narrow, every small bump becomes a mode and you get too many micro-clusters.

The heuristic bandwidth estimation following Silverman’s rule is often used as a starting point, then refined empirically.

Python Implementation with scikit-learn

Basic example with MeanShift

Here is a complete implementation using scikit-learn:

from sklearn.cluster import MeanShift, estimate_bandwidth
from sklearn.datasets import make_blobs
import matplotlib.pyplot as plt
import numpy as np

# Generate synthetic data
X, y_true = make_blobs(
    n_samples=500,
    centers=4,
    cluster_std=1.0,
    random_state=42
)

# Automatic bandwidth estimation
bandwidth = estimate_bandwidth(X, quantile=0.2, n_samples=200)
print(f"Estimated bandwidth: {bandwidth:.3f}")

# Run Mean Shift
ms = MeanShift(bandwidth=bandwidth, bin_seeding=True)
ms.fit(X)
labels = ms.labels_
cluster_centers = ms.cluster_centers_

n_clusters = len(set(labels))
print(f"Number of clusters detected: {n_clusters}")

# Visualization
plt.figure(figsize=(10, 6))
colors = plt.cm.tab10(np.linspace(0, 1, n_clusters))

for i, color in enumerate(colors):
    mask = labels == i
    plt.scatter(X[mask, 0], X[mask, 1], c=[color],
                label=f"Cluster {i}", alpha=0.6, s=30)

# Draw centers
plt.scatter(cluster_centers[:, 0], cluster_centers[:, 1],
            c="red", marker="X", s=200, edgecolors="black",
            label="Centers", zorder=5)

plt.title(f"Mean Shift — {n_clusters} clusters detected (bandwidth={bandwidth:.2f})")
plt.legend()
plt.tight_layout()
plt.show()

Comparison with K-Means

Here’s how to compare Mean Shift and K-Means on the same dataset:

from sklearn.cluster import KMeans

# K-Means (assuming we know the right K)
kmeans = KMeans(n_clusters=4, random_state=42, n_init="auto")
kmeans_labels = kmeans.fit_predict(X)

# Mean Shift
ms = MeanShift(bandwidth=bandwidth, bin_seeding=True)
ms_labels = ms.fit_predict(X)

# Visual comparison
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# K-Means
ax1 = axes[0]
for i in range(4):
    mask = kmeans_labels == i
    ax1.scatter(X[mask, 0], X[mask, 1], alpha=0.6, s=30,
                label=f"Cluster {i}")
ax1.scatter(kmeans.cluster_centers_[:, 0], kmeans.cluster_centers_[:, 1],
            c="red", marker="X", s=200, edgecolors="black")
ax1.set_title("K-Means (K=4 fixed)")
ax1.legend()

# Mean Shift
ax2 = axes[1]
unique_labels = set(ms_labels)
for i, label in enumerate(unique_labels):
    mask = ms_labels == label
    ax2.scatter(X[mask, 0], X[mask, 1], alpha=0.6, s=30,
                label=f"Cluster {i}")
ax2.scatter(ms.cluster_centers_[:, 0], ms.cluster_centers_[:, 1],
            c="red", marker="X", s=200, edgecolors="black")
ax2.set_title(f"Mean Shift ({len(unique_labels)} automatic clusters)")
ax2.legend()

plt.tight_layout()
plt.show()

Convergence visualization

We can observe the convergence path of each point toward its mode:

# Run Mean Shift with center computation
ms = MeanShift(bandwidth=bandwidth, bin_seeding=True)
ms.fit(X)

# Visualize convergence for a subsample
sample_indices = np.random.choice(len(X), size=80, replace=False)

fig, ax = plt.subplots(figsize=(8, 6))

# Data points
ax.scatter(X[:, 0], X[:, 1], c=ms.labels_, cmap="tab10",
           alpha=0.3, s=20)

# Convergence lines
for idx in sample_indices:
    point = X[idx]
    center = ms.cluster_centers_[ms.labels_[idx]]
    ax.plot([point[0], center[0]], [point[1], center[1]],
            color="gray", alpha=0.15, linewidth=0.8)

# Final centers
ax.scatter(ms.cluster_centers_[:, 0], ms.cluster_centers_[:, 1],
           c="red", marker="X", s=200, edgecolors="black", zorder=5)
ax.set_title("Convergence paths toward modes")
plt.tight_layout()
plt.show()

Mean Shift Hyperparameters

Hyperparameter selection is crucial to get good results with Mean Shift.

bandwidth

This is the most important parameter. It controls the radius of the Parzen window used to estimate the density.

• Too small: every small group of points becomes its own cluster → over-segmentation
• Too large: all points merge into a single cluster → under-segmentation
• Ideal value: corresponds to the characteristic scale of the structures you want to detect

The heuristic estimate_bandwidth(X, quantile=0.2) computes the quantile of pairwise distances, providing a good starting point. Then adjust the quantile (between 0.1 and 0.5) based on the desired granularity.

bin_seeding

When bin_seeding=True, the algorithm discretizes the space into bins and only initializes clustering seeds once per occupied bin. This significantly reduces computation time, especially on large datasets. The bin size is equal to the bandwidth.

min_bin_freq

This parameter controls the minimum number of points a bin must contain to be considered a valid seed. With min_bin_freq=1 (default value), every occupied bin becomes a seed. Increasing this value ignores less dense regions and speeds up computation, at the cost of potentially missing small clusters.

cluster_all

Determines the treatment of orphan points that do not clearly converge toward a mode.

• cluster_all=True (default): all points are assigned, including borderline cases that are attached to the nearest cluster
• cluster_all=False: unassigned points receive the label -1

n_jobs

Number of processors used in parallel. On a modern processor, n_jobs=-1 uses all available cores and significantly reduces the density estimation computation time.

Advantages and Limitations of Mean Shift

Advantages

1. No need to fix the number of clusters — the number of modes is determined automatically by the data structure.
2. Clusters of arbitrary shape — basins of attraction are not constrained to be spherical or convex.
3. Robust to outliers — isolated points converge to their own mode or are identified as orphans (depending on cluster_all).
4. Solid theoretical foundation — based on kernel density estimation, with proven convergence guarantees.
5. Interpretability — each cluster corresponds to an identifiable mode of the underlying distribution.

Limitations

1. High computational cost — naive complexity is O(n²·d·T) where T is the number of iterations. For n > 10,000 points, this becomes prohibitive without bin_seeding.
2. Bandwidth sensitivity — choosing an inappropriate bandwidth gives poor results. There is no universal automatic method for selecting it.
3. Difficulty in high dimensions — in high-dimensional space, kernel density estimation suffers from the curse of dimensionality: all points become equally close and the notion of local density loses its meaning.
4. Clusters of very different densities — a single bandwidth may perform poorly if clusters have very heterogeneous densities (one very dense cluster and another very sparse one).
5. No probabilistic model — unlike Gaussian Mixture Models, Mean Shift does not provide membership probabilities or likelihood.

4 Concrete Use Cases of Mean Shift

Case 1 — Image segmentation

Mean Shift is historically the reference algorithm for image segmentation. Each pixel is represented as a point in a space combining spatial coordinates (x, y) and color (R, G, B). The bandwidth controls the trade-off between color fidelity and spatial regularity.

from sklearn.cluster import MeanShift
from sklearn.cluster import estimate_bandwidth
import numpy as np
from skimage import data

# Load an image
image = data.astronaut()
height, width, channels = image.shape

# Prepare data: concatenate position and color
X_spatial = np.mgrid[0:height, 0:width].reshape(2, -1).T
X_color = image.reshape(-1, 3).astype(float)
X = np.hstack([X_spatial * 0.3, X_color])

# Mean Shift
bandwidth = estimate_bandwidth(X, quantile=0.05, n_samples=500)
ms = MeanShift(bandwidth=bandwidth, bin_seeding=True)
ms.fit(X[::10])  # sample to speed up

# Reconstruct the segmented image
labels_full = ms.predict(X)
segmented = labels_full.reshape(height, width)
print(f"Image segmented into {len(set(labels_full))} regions")

Case 2 — GPS trajectory analysis

In geolocation, Mean Shift is used to identify frequently visited places (points of interest) from GPS traces. Points naturally cluster around visited locations, forming spatial density modes.

# Simulated GPS coordinates
gps_data = np.vstack([
    np.random.randn(200, 2) * 0.001 + [48.8566, 2.3522],
    np.random.randn(150, 2) * 0.002 + [48.8738, 2.2950],
    np.random.randn(100, 2) * 0.001 + [48.8448, 2.3736],
    np.random.randn(80, 2) * 0.003 + [48.8600, 2.3400],
])

bandwidth = estimate_bandwidth(gps_data, quantile=0.15)
ms = MeanShift(bandwidth=bandwidth, bin_seeding=True)
ms.fit(gps_data)

for i, center in enumerate(ms.cluster_centers_):
    print(f"Frequent location #{i+1}: lat={center[0]:.4f}, lon={center[1]:.4f} "
          f"({np.sum(ms.labels_==i)} points)")

Case 3 — Community detection in social networks

By representing users as points in a latent space (embeddings from Node2Vec, for example), Mean Shift automatically detects communities without imposing their number. Groups of users with similar profiles form dense modes.

from sklearn.decomposition import PCA

# X_embedded: embedding matrix (n_users, 300)
pca = PCA(n_components=50)
X_reduced = pca.fit_transform(X_embedded)

bandwidth = estimate_bandwidth(X_reduced, quantile=0.1, n_samples=500)
ms = MeanShift(bandwidth=bandwidth, bin_seeding=True, n_jobs=-1)
ms.fit(X_reduced)

n_communities = len(set(ms.labels_))
print(f"{n_communities} communities detected")
for c in range(n_communities):
    print(f"  Community {c}: {np.sum(ms.labels_==c)} users")

Case 4 — Climate regime identification

In climatology, Mean Shift is used to identify atmospheric circulation regimes from geopotential data. Each mode corresponds to a recurring circulation pattern, which helps understand natural climate variability.

# Simulated climate data
X_climate = np.vstack([
    np.random.randn(500, 3) * 0.5 + [1, 0.5, -0.3],
    np.random.randn(350, 3) * 0.4 + [-1, -0.5, 0.3],
    np.random.randn(250, 3) * 0.6 + [0.5, -1, 0.8],
    np.random.randn(400, 3) * 0.3 + [-0.3, 1, -0.5],
])

bandwidth = estimate_bandwidth(X_climate, quantile=0.2)
ms = MeanShift(bandwidth=bandwidth, bin_seeding=True)
ms.fit(X_climate)

print(f"{len(set(ms.labels_))} climate regimes identified")
for i, center in enumerate(ms.cluster_centers_):
    n_points = np.sum(ms.labels_ == i)
    print(f"  Regime {i+1}: center={center.round(2)}, {n_points} days")

When to choose Mean Shift over another algorithm?

Choose Mean Shift when:

• The number of clusters is unknown and you want the algorithm to determine it itself
• You are working with data of moderate dimensionality (d < 50)
• Clusters potentially have non-spherical shapes
• The dataset is of reasonable size (n < 10,000 without bin_seeding, or n < 100,000 with bin_seeding)
• You need a density-based method but DBSCAN fails to capture the structure

Avoid Mean Shift when:

• Data is in very high dimensions (> 100) — prefer preliminary dimensionality reduction
• The dataset exceeds a few hundred thousand points — the computational cost becomes prohibitive
• Clusters have very heterogeneous densities — DBSCAN or HDBSCAN will be more suitable
• You need probabilistic assignments — turn to Gaussian Mixture Models

See Also

• Master Maximal Arrangements in Python: Complete Guide and Tips for Developers
• Finding the Median of a Data Stream in Python — Interview Question