Détectez des zones denses avec DBSCAN

Sur une carte de Rennes couverte de photos géolocalisées, certains endroits concentrent des dizaines de clichés (une place, un monument), d’autres n’en ont qu’un, pris au hasard d’une promenade. Combien y a-t-il de lieux d’intérêt ? Personne ne le sait. Et que faire des photos isolées ? Les ranger de force dans le groupe le plus proche fausserait tout.

k-means, vu au chapitre précédent, est mal armé pour ce problème. Ce chapitre présente DBSCAN, un algorithme qui regroupe les points par densité, trouve seul le nombre de groupes et sait dire « ce point n’appartient à aucun groupe ». Notre fil rouge : dix photos placées sur un plan, en centaines de mètres, que vous pourrez traiter à la main, avant de passer aux vraies coordonnées GPS.

Pourquoi k-means ne suffit pas

k-means a trois défauts face à des données spatiales réelles :

il faut lui donner le nombre de groupes $k$ , alors qu’on cherche justement combien de lieux d’intérêt il y a ;
il produit des groupes convexes, en boule autour d’un centre ; un groupe de photos qui longe une rivière, ou qui forme un anneau autour d’une place, lui échappe ;
il affecte chaque point à un groupe : la photo isolée prise dans un jardin privé est rattachée au monument le plus proche, dont elle déplace le centre.

k-means a aussi des qualités : il est très rapide (son coût est linéaire en nombre de points) et n’a pas besoin de calculer toutes les distances entre paires de points. Mais pour trouver des zones denses dans un nuage bruité, il faut changer d’idée.

L’idée : un groupe est une zone dense

DBSCAN (density-based spatial clustering of applications with noise) définit un groupe comme une région où les points sont serrés, entourée de régions clairsemées. Il repose sur deux paramètres :

eps ( $\varepsilon$ ) : le rayon du voisinage d’un point ;
min_samples : le nombre minimal de points dans ce voisinage (le point lui-même compris) pour que la zone soit dite dense.

On classe alors chaque point dans l’une de trois catégories :

un point cœur (core point) a au moins min_samples points dans son $\varepsilon$ -voisinage ;
un point de bordure (border point) n’est pas cœur, mais se trouve dans le voisinage d’un point cœur ;
un point de bruit (noise) n’est ni l’un ni l’autre.

Un groupe est formé de points cœurs connectés de proche en proche (chacun dans le voisinage d’un autre), plus les points de bordure qui leur sont rattachés. L’algorithme est simple : on prend un point cœur non encore visité, on ouvre un nouveau groupe, on y ajoute tous ses voisins, et l’on étend récursivement à partir de chaque voisin qui est lui-même cœur. On recommence jusqu’à ce qu’il ne reste plus de point cœur non visité. Ce qui reste hors de tout groupe est du bruit, étiqueté −1 dans scikit-learn.

Dix photos ont pour positions (en centaines de mètres) : A (1 ; 1), B (2 ; 1), C (1 ; 2), D (2 ; 2), E (3,2 ; 2), F (6 ; 5), G (7 ; 5), H (6,5 ; 6), I (4,5 ; 5,5) et J (9 ; 1). On prend eps = 1,5 et min_samples = 3.

Point	Voisins à 1,5 ou moins (lui compris)	Nombre	Statut
A, B, C	A, B, C, D	4	cœur
D	A, B, C, D, E	5	cœur
E	D, E	2	bordure (voisin du cœur D)
F, G, H	F, G, H	3	cœur
I	I	1	bruit
J	J	1	bruit

On obtient deux groupes : $\lbrace A, B, C, D, E \rbrace$ et $\lbrace F, G, H \rbrace$ , et deux points de bruit, I et J. Le point I, à 1,58 de F, est juste hors de portée.

Les cercles en pointillés montrent le voisinage de rayon eps autour de D et de F. Avec un k-means à deux groupes, I et J auraient été rattachés de force à l’un des groupes, et J aurait tiré le centre de son groupe vers la droite.

Réglez eps et min_samples

DBSCAN n’a pas besoin du nombre de groupes, mais ses deux paramètres comptent beaucoup.

min_samples fixe ce qu’on appelle « dense ». Plus il est grand, plus les groupes doivent être fournis, et plus il y a de bruit. Une règle courante est de le prendre au moins égal à la dimension plus un, souvent deux fois la dimension pour des données bruitées, puis d’ajuster selon le sens métier : pour un lieu d’intérêt, combien de photos faut-il au minimum ?

eps fixe l’échelle. Trop petit, presque tout devient du bruit ; trop grand, des groupes distincts fusionnent. La méthode classique est le graphique des k-distances : pour chaque point, on calcule la distance à son $k$ -ième plus proche voisin (avec $k$ = min_samples − 1), on trie ces distances et on les trace. La courbe reste plate pour les points des zones denses, puis se redresse brusquement pour les points isolés : on choisit eps au niveau de ce coude.

La courbe reste plate pour les points des deux groupes, puis se redresse pour E (bordure), I et J (bruit). Tout eps pris dans la bande verte, comme 1,5, donne le résultat du tableau.

import numpy as np
from sklearn.cluster import DBSCAN
from sklearn.neighbors import NearestNeighbors

X = np.array([[1, 1], [2, 1], [1, 2], [2, 2], [3.2, 2],
              [6, 5], [7, 5], [6.5, 6], [4.5, 5.5], [9, 1]])

# Graphique des k-distances : distance au 2e voisin (le point lui-même est le voisin 0)
dist, _ = NearestNeighbors(n_neighbors=3).fit(X).kneighbors(X)
print(np.sort(dist[:, 2]).round(2))

db = DBSCAN(eps=1.5, min_samples=3).fit(X)
print(db.labels_)                    # [ 0  0  0  0  0  1  1  1 -1 -1] : -1 = bruit
print(db.core_sample_indices_)       # indices des points cœurs

Attention aux unités : les coordonnées GPS

Comme k-means, DBSCAN repose sur une distance, et l’unité de cette distance détermine le sens de eps. C’est particulièrement piégeux avec des coordonnées GPS.

Une latitude et une longitude s’expriment en degrés, et un degré ne mesure pas la même longueur dans les deux directions. Un degré de latitude vaut environ 111 km partout ; un degré de longitude vaut 111 km multipliés par le cosinus de la latitude, soit environ 74 km à Rennes (latitude 48°). Un eps de 0,0003 degré correspond donc à environ 33 m dans le sens nord-sud, mais seulement 22 m dans le sens est-ouest : le voisinage n’est pas un cercle, mais une ellipse.

Un eps en degrés donne un voisinage aplati d'est en ouest à la latitude de Rennes ; la distance haversine (ou une projection en mètres) redonne un vrai cercle.

Pour bien faire, on a deux options :

projeter les coordonnées dans un système métrique (par exemple Lambert 93 en France), puis utiliser la distance euclidienne en mètres ;
utiliser la distance haversine (la distance sur la sphère terrestre), qui attend des coordonnées en radians et renvoie une distance en radians, qu’on convertit avec le rayon de la Terre.

import numpy as np
from sklearn.cluster import DBSCAN

R_TERRE = 6_371_000                               # rayon moyen de la Terre, en mètres
coords = np.radians(photos[["lat", "long"]].to_numpy())
db = DBSCAN(eps=30 / R_TERRE,                     # un voisinage de 30 mètres
            min_samples=7, metric="haversine", algorithm="ball_tree").fit(coords)

À l’échelle d’une ville et pour une première analyse, la distance euclidienne sur les degrés reste une approximation acceptable, à condition de savoir ce qu’elle mesure. C’est ce que nous ferons en TP (chapitre 10).

Forces et limites de DBSCAN

	k-means	CAH (Ward)	DBSCAN
Nombre de groupes	à fixer ( $k$ )	choisi en coupant le dendrogramme	trouvé par l’algorithme
Forme des groupes	convexe	plutôt compacte	quelconque
Points isolés	rattachés à un groupe	rattachés à un groupe	étiquetés comme bruit
Paramètres	$k$ , initialisation	critère de fusion	`eps`, `min_samples`
Coût	faible (linéaire)	élevé (au moins quadratique)	modéré avec un index spatial
Point faible	groupes non convexes, bruit	grands volumes	densités variables

Données simulées. k-means découpe l'espace en deux zones convexes autour de ses centres, quitte à couper chaque croissant ; DBSCAN suit les zones denses, quelle que soit leur forme, et isole le bruit.

La limite principale de DBSCAN est qu’il n’utilise qu’une seule densité pour tout le jeu de données. Or un centre-ville, où les photos se pressent, et un grand parc, où elles sont plus espacées, demanderaient des eps différents : un eps adapté au centre émiette le parc en bruit, un eps adapté au parc fusionne tout le centre en un seul groupe. Des variantes comme OPTICS et HDBSCAN traitent ce problème en explorant toute une gamme de densités.

Données simulées. Aucun eps unique ne convient : petit, il sépare les lieux du centre mais émiette le parc ; grand, il trouve le parc mais fond le centre en un seul groupe.

Côté coût, DBSCAN doit trouver les voisins de chaque point. Avec un index spatial (arbre k-d ou ball tree), cela se fait en $O(n \log n)$ en moyenne ; sans index, le coût devient quadratique. En très grande dimension, le fléau de la dimension (chapitre 4) rend la notion même de densité fragile.

En résumé

k-means impose le nombre de groupes, des formes convexes et l’affectation de tous les points, ce qui convient mal aux données spatiales bruitées.
DBSCAN regroupe par densité avec deux paramètres : le rayon eps et le seuil min_samples (point lui-même compris).
Un point est cœur, de bordure ou bruit ; un groupe est un ensemble de points cœurs connectés et de leurs bordures ; le bruit est étiqueté −1.
On règle eps avec le graphique des k-distances, et min_samples selon la dimension et le sens métier.
Sur des coordonnées GPS, un eps en degrés n’a pas la même longueur en latitude et en longitude : on projette ou l’on utilise la distance haversine.
DBSCAN ne gère qu’une densité à la fois ; OPTICS et HDBSCAN lèvent cette limite.

Exercices

Exercice · Classer des points à la main

facile

Sur une droite, sept points sont placés en 0 ; 0,5 ; 1 ; 1,4 ; 3 ; 6 et 6,4. On applique DBSCAN avec eps = 0,6 et min_samples = 3. Donnez le statut de chaque point et les groupes obtenus.

Voir la solution

On liste les voisins à distance 0,6 ou moins, le point compris :

0 : $\lbrace 0 ;\, 0{,}5 \rbrace$ , 2 points, pas cœur ;
0,5 : $\lbrace 0 ;\, 0{,}5 ;\, 1 \rbrace$ , 3 points, cœur ;
1 : $\lbrace 0{,}5 ;\, 1 ;\, 1{,}4 \rbrace$ , 3 points, cœur ;
1,4 : $\lbrace 1 ;\, 1{,}4 \rbrace$ , 2 points, pas cœur ;
3 : seul, pas cœur ;
6 et 6,4 : $\lbrace 6 ;\, 6{,}4 \rbrace$ , 2 points chacun, pas cœurs.

Les cœurs 0,5 et 1 sont voisins : ils forment un groupe, auquel s’ajoutent les points de bordure 0 (voisin de 0,5) et 1,4 (voisin de 1). Le groupe est $\lbrace 0 ;\, 0{,}5 ;\, 1 ;\, 1{,}4 \rbrace$ . Les points 3, 6 et 6,4 sont du bruit : 6 et 6,4 sont proches l’un de l’autre, mais deux points ne suffisent pas à former une zone dense avec min_samples = 3.

En code

Python
R

import numpy as np
from sklearn.cluster import DBSCAN

x = np.array([0, 0.5, 1, 1.4, 3, 6, 6.4]).reshape(-1, 1)
db = DBSCAN(eps=0.6, min_samples=3).fit(x)

# Taille de chaque voisinage (le point compris)
print((np.abs(x - x.T) <= 0.6).sum(axis=1))     # [2 3 3 2 1 2 2]
print(db.labels_)                                # [ 0  0  0  0 -1 -1 -1]
print(x[db.core_sample_indices_].ravel())        # [0.5 1. ] : les deux points cœurs

# DBSCAN minimal (base R) : -1 = bruit, groupes numérotés à partir de 0
dbscan_simple <- function(X, eps, min_samples) {
  voisins <- as.matrix(dist(X)) <= eps          # voisinage, le point compris
  coeur <- rowSums(voisins) >= min_samples
  etiquette <- rep(-1, nrow(voisins)); groupe <- -1
  for (i in which(coeur)) if (etiquette[i] == -1) {
    groupe <- groupe + 1; a_voir <- i
    while (length(a_voir) > 0) {
      j <- a_voir[1]; a_voir <- a_voir[-1]
      if (etiquette[j] != -1) next
      etiquette[j] <- groupe
      if (coeur[j]) a_voir <- c(a_voir, which(voisins[j, ] & etiquette == -1))
    }
  }
  list(etiquettes = etiquette, coeur = coeur)
}

x <- c(0, 0.5, 1, 1.4, 3, 6, 6.4)
res <- dbscan_simple(x, eps = 0.6, min_samples = 3)
rowSums(as.matrix(dist(x)) <= 0.6)   # 2 3 3 2 1 2 2 : taille des voisinages
res$etiquettes                       # 0 0 0 0 -1 -1 -1
x[res$coeur]                         # 0.5 1 : les deux points cœurs

Exercice · Faire varier les paramètres

moyen

Reprenez les dix photos du chapitre (A à J). Que devient le résultat avec eps = 1,5 et min_samples = 4 ? Et avec eps = 1,0 et min_samples = 3 ? Commentez.

Voir la solution

eps = 1,5 et min_samples = 4. A, B et C ont 4 voisins et D en a 5 : ils restent cœurs. F, G et H n’ont que 3 voisins : ils ne sont plus cœurs, et aucun cœur ne les atteint. Résultat : un seul groupe $\lbrace A, B, C, D, E \rbrace$ (E reste en bordure de D), et cinq points de bruit : F, G, H, I et J. Le petit groupe de trois photos a disparu.

eps = 1,0 et min_samples = 3. Les distances de 1 exactement comptent (A-B, A-C, B-D, C-D), mais pas la diagonale A-D ( $\sqrt 2$ ) ni D-E (1,2). A, B, C et D ont chacun 3 voisins : ils sont cœurs. F et G sont à distance 1, mais H est à 1,12 de chacun : F et G n’ont que 2 voisins. Résultat : un seul groupe $\lbrace A, B, C, D \rbrace$ , et six points de bruit (E, F, G, H, I, J).

Dans les deux cas, durcir un seul paramètre fait basculer tout un groupe dans le bruit : DBSCAN est sensible à ses réglages, d’où l’intérêt du graphique des k-distances et d’un seuil min_samples justifié par le métier.

En code

Python
R

import numpy as np
from sklearn.cluster import DBSCAN

noms = np.array(list("ABCDEFGHIJ"))
X = np.array([[1, 1], [2, 1], [1, 2], [2, 2], [3.2, 2],
              [6, 5], [7, 5], [6.5, 6], [4.5, 5.5], [9, 1]])

for eps, min_samples in [(1.5, 3), (1.5, 4), (1.0, 3)]:
    db = DBSCAN(eps=eps, min_samples=min_samples).fit(X)
    print(eps, min_samples, db.labels_, "bruit :", "".join(noms[db.labels_ == -1]))
# 1.5 3 [ 0  0  0  0  0  1  1  1 -1 -1] bruit : IJ
# 1.5 4 [ 0  0  0  0  0 -1 -1 -1 -1 -1] bruit : FGHIJ
# 1.0 3 [ 0  0  0  0 -1 -1 -1 -1 -1 -1] bruit : EFGHIJ

# DBSCAN minimal (base R) : -1 = bruit, groupes numérotés à partir de 0
dbscan_simple <- function(X, eps, min_samples) {
  voisins <- as.matrix(dist(X)) <= eps          # voisinage, le point compris
  coeur <- rowSums(voisins) >= min_samples
  etiquette <- rep(-1, nrow(voisins)); groupe <- -1
  for (i in which(coeur)) if (etiquette[i] == -1) {
    groupe <- groupe + 1; a_voir <- i
    while (length(a_voir) > 0) {
      j <- a_voir[1]; a_voir <- a_voir[-1]
      if (etiquette[j] != -1) next
      etiquette[j] <- groupe
      if (coeur[j]) a_voir <- c(a_voir, which(voisins[j, ] & etiquette == -1))
    }
  }
  list(etiquettes = etiquette, coeur = coeur)
}

X <- matrix(c(1, 1,  2, 1,  1, 2,  2, 2,  3.2, 2,
              6, 5,  7, 5,  6.5, 6,  4.5, 5.5,  9, 1), ncol = 2, byrow = TRUE,
            dimnames = list(LETTERS[1:10], NULL))
for (reglage in list(c(1.5, 3), c(1.5, 4), c(1.0, 3))) {
  e <- dbscan_simple(X, eps = reglage[1], min_samples = reglage[2])$etiquettes
  cat(reglage, ":", e, "| bruit :", rownames(X)[e == -1], "\n")
}
# 1.5 3 : 0 0 0 0 0 1 1 1 -1 -1 | bruit : I J
# 1.5 4 : 0 0 0 0 0 -1 -1 -1 -1 -1 | bruit : F G H I J
# 1 3 : 0 0 0 0 -1 -1 -1 -1 -1 -1 | bruit : E F G H I J

Exercice · Convertir un rayon en degrés

moyen

Vous voulez un voisinage d’environ 50 m autour de chaque photo prise à Rennes (latitude 48,1°), en travaillant directement sur les degrés. Quel eps choisir pour que 50 m soient respectés dans le sens nord-sud ? Que mesure alors ce même eps dans le sens est-ouest ? Quel eps passeriez-vous à DBSCAN avec la métrique haversine ?

Voir la solution

Un degré de latitude mesure environ $111{,}32$ km. Pour 50 m : $\text{eps} = 0{,}05 / 111{,}32 \approx 0{,}00045$ degré.

Dans le sens est-ouest, un degré de longitude mesure $111{,}32 \times \cos(48{,}1^ rc) \approx 111{,}32 \times 0{,}668 \approx 74{,}3$ km. Le même eps de 0,00045 degré y représente donc environ $0{,}00045 \times 74\,300 \approx 33$ m seulement. Le voisinage est une ellipse, plus étroite d’est en ouest.

Avec la métrique haversine, les coordonnées sont en radians et eps est un angle : $\text{eps} = 50 / 6\,371\,000 \approx 7{,}85 \times 10^{-6}$ radian. Le voisinage est alors un vrai disque de 50 m dans toutes les directions.

En code

Python
R

import numpy as np

km_par_degre = 111.32                      # un degré de latitude
lat = np.radians(48.1)                     # Rennes
rayon_m = 50

eps_deg = rayon_m / 1000 / km_par_degre
print(round(eps_deg, 5))                                  # 0.00045 degré
km_longitude = km_par_degre * np.cos(lat)
print(round(km_longitude, 1), round(eps_deg * km_longitude * 1000))   # 74.3 km ; 33 m
print(f"{rayon_m / 6_371_000:.3g}")                       # 7.85e-06 radian (haversine)

km_par_degre <- 111.32                     # un degré de latitude
lat <- 48.1 * pi / 180                     # Rennes, en radians
rayon_m <- 50

eps_deg <- rayon_m / 1000 / km_par_degre
round(eps_deg, 5)                                        # 0.00045 degré
km_longitude <- km_par_degre * cos(lat)
round(c(km_longitude, eps_deg * km_longitude * 1000), 1) # 74.3 km ; 33.4 m
signif(rayon_m / 6371000, 3)                             # 7.85e-06 radian (haversine)