Détection de ruptures

En analyse statistique, le problème de détection de ruptures (ou détection de points de changement) est un problème de régression ayant pour but d'estimer les instants où un signal présente des changements dans la distribution. Ces instants sont matérialisés sur les deux figures par des lignes verticales bleues. Classiquement, on réalise de la détection de rupture pour un signal ayant des changements dans la moyenne. De manière plus générale, on réalise de la détection de ruptures pour un signal ayant des changements dans la distribution (par exemple, dans la moyenne et la variance).

La détection de ruptures peut s'appliquer à un signal sonore d'une émission dont on souhaite estimer les instants où l'on change de situations, à la détection d'attaques informatiques (variations de flux réseaux)^[1], ou encore au contrôle qualité.

Cet article traite du problème de détection de ruptures rétrospective (dite offline) où l'on dispose de l'ensemble des données du signal. Ce contexte est différent d'une détection temps réel (dite online) où les données du signal arrivent à la volée, mais moins à même de détecter précisément l'instant de rupture.

Formalisation du problème

Soit ${X_{t}}_{t \in [[1, T]]}$ un signal réel, provenant d'observations recueillies au cours des instants $1, \dots, T$ , présentant des changements dans la distribution. En notant $P_{X_{i}}$ la loi de probabilité de $X_{i}$ , la distribution de ${X_{t}}_{t \in [[1, T]]}$ vérifie :

\begin{matrix} P_{X_{τ_{0}^{*}}} = \dots = P_{X_{τ_{1}^{*} - 1}} \neq P_{X_{τ_{1}^{*}}} = \dots = P_{X_{τ_{2}^{*} - 1}} \neq \dots \neq P_{X_{τ_{D^{*} - 1}^{*}}} = \dots = P_{X_{τ_{D^{*}}^{*} - 1}}, \end{matrix}

avec $τ_{0}^{*}, \dots, τ_{D^{*}}^{*}$ étant les vrais instants de ruptures (on note $m^{*} = {τ_{0}^{*}, \dots, τ_{D^{*}}^{*}}$ ( $D^{*}$ est le vrai nombre de segments) avec la convention $τ_{0}^{*} = 1$ et $τ_{D^{*}}^{*} = T + 1$ ). On cherche à estimer ces instants de ruptures à l'aide d'un algorithme.

Dans le cas de la détection de rupture dans la moyenne, le modèle est :

\begin{matrix} \forall t \in [[1, T]], X_{t} = s_{t} + ϵ_{t}, \end{matrix}

avec ${s_{t}}_{t \in [[1, T]]}$ est la fonction de régression et ${ϵ_{t}}_{t \in [[1, T]]}$ est un bruit d'espérance nulle et de variance $σ^{2}$ . La fonction de régression ${s_{t}}_{t \in [[1, T]]}$ est supposée constante par morceaux avec des discontinuités à chaque vrai instant de ruptures ${τ_{i}^{*}}_{i \in [[0, D^{*}]]}$ .

Dans le cas de la détection de ruptures dans la distribution, on recode les observations initiales ${X_{t}}_{t \in [[1, T]]}$ par de nouvelles observations ${Y_{t}}_{t \in [[1, T]]}$ définies par $Y_{t} = k_{X_{t}} = k (X_{t}, .)$ où $k$ est un noyau symétrique et semi-défini positif (par exemple $\forall x, y \in ℝ, k^{1} (x, y) = x * y$ : $k^{1}$ est le noyau linéaire ; autre exemple $\forall x, y \in ℝ, k_{δ}^{2} (x, y) = \exp (- \frac{(x - y)^{2}}{2 δ^{2}})$ : $k_{δ}^{2}$ est le noyau Gaussien de paramètre $δ > 0$ ). Pour un noyau symétrique et semi-défini positif $k$ , le théorème de Moore-Aronszahn assure l'existence d'un espace de Hilbert à noyau reproduisant $ℋ$ de noyau reproduisant $k$ .

Le modèle est :

\begin{matrix} \forall t \in [[1, T]], Y_{t} = μ_{t}^{*} + ϵ_{t}, \end{matrix}

avec ${μ_{t}^{*}}_{t \in [[1, T]]}$ est la fonction de régression et ${ϵ_{t}}_{t \in [[1, T]]}$ est un bruit d'espérance nulle et de variance $σ^{2}$ . De plus, $\forall t \in [[1, T]], Y_{t}, μ_{t}^{*}, ϵ_{t}$ appartiennent à $ℋ$ . La fonction de régression ${μ_{t}^{*}}_{t \in [[1, T]]}$ est supposée constante par morceaux avec des discontinuités à chaque vrai instant de ruptures ${τ_{i}^{*}}_{i \in [[0, D^{*}]]}$ .

Méthodes existantes

Le problème de détection de ruptures peut être vu comme un problème de sélection de modèle^[2]: chaque segmentation candidate (i.e. liste d'instants de ruptures candidats) est reliée à un modèle qu'il faut choisir. Nous présentons deux approches : l'une utilisant le principe de la programmation dynamique et l'autre l'heuristique de segmentation binaire dans le cadre classique de la détection de ruptures dans la moyenne puis dans le cas général de la détection de ruptures dans la distribution. Pour chacune de ces deux méthodes, nous présentons les algorithmes permettant de calculer une segmentation en $D$ segments ${\hat{m}}_{D}$ optimisant un critère (ici, le risque empirique).

Notations

Nous présentons les notations communes aux deux méthodes puis celles qui leur sont spécifiques :

Notation commune

$ℳ (D)$ l'ensemble des segmentions $m = {{τ_{0}, \dots, τ_{1} - 1}, {τ_{1}, \dots, τ_{2} - 1}, \dots, {τ_{D - 1}, \dots, τ_{D} - 1}}$ en $D$ segments $(avec la convention τ_{0} = 1, τ_{D} = T + 1)$ .

Notations pour le cas de la détection de ruptures dans la moyenne

$X = (X_{1}, \dots, X_{T})$ et $s = (s_{1}, \dots, s_{T})$ .
$S_{m}^{ℝ} = {u \in ℝ^{T} : u_{τ_{0}} = \dots = u_{τ_{1} - 1} \neq u_{τ_{1}} = \dots = u_{τ_{2} - 1} \neq \dots \neq u_{τ_{D - 1}} = \dots u_{τ_{D} - 1}}$ : $S_{m}^{ℝ}$ est l'ensemble des vecteurs $u \in ℝ^{T}$ constants par morceaux sur les segments de $m$ .
Estimateur du risque empirique : ${\hat{s}}_{m} = \underset{u \in S_{m}^{ℝ}}{a r g m i n} ‖ X - u ‖_{2}^{2}$ .

Notations pour le cas de la détection de ruptures dans la distribution

$Y = (Y_{1}, \dots, Y_{T})$ et $μ^{*} = (μ_{1}^{*}, \dots, μ_{T}^{*})$ .
$Y_{s}^{e} = {Y_{s}, \dots, Y_{e}}$ pour $s < e$ avec $(s, e) \in [[1, T]]^{2}$ .
$S_{m}^{ℋ} = {u \in ℋ^{T} : u_{τ_{0}} = \dots = u_{τ_{1} - 1} \neq u_{τ_{1}} = \dots = u_{τ_{2} - 1} \neq \dots \neq u_{τ_{D - 1}} = \dots u_{τ_{D} - 1}}$ : $S_{m}^{ℋ}$ est l'ensemble des vecteurs $u \in ℋ^{T}$ constants par morceaux sur les segments de $m$ .
Norme dans $ℋ^{T}$ : $\forall h \in ℋ^{T}, ‖ h ‖_{ℋ^{T}}^{2} = \sum_{i = 1}^{T} ‖ h_{i} ‖_{ℋ}^{2}$ .
Estimateur du risque empirique : ${\hat{μ}}_{m} = \underset{u \in S_{m}^{ℋ}}{a r g m i n} ‖ Y - u ‖_{ℋ^{T}}^{2}$ .

Les méthodes proposées ci-dessous utilisent le risque empirique comme critère à minimiser ( $‖ X - {\hat{s}}_{m} ‖_{2}^{2}$ pour la détection de ruptures dans la moyenne ; $‖ Y - {\hat{μ}}_{m} ‖_{ℋ^{T}}^{2}$ pour la détection de ruptures dans la distribution). Pour le noyau linéaire $k^{1}$ , $‖ X - {\hat{s}}_{m} ‖_{2}^{2} = ‖ Y - {\hat{μ}}_{m} ‖_{ℋ^{T}}^{2}$ , les méthodes utilisées dans le cadre classique se déduisent de celles du cas des noyaux par le biais de $k^{1}$ : donc on ne présentera les algorithmes que dans le cas des noyaux.

Programmation dynamique

La méthode de la programmation dynamique utilise le principe d'optimalité de Bellman : toute solution optimale s'appuie elle-même sur des sous-problèmes résolus localement de façon optimale. On utilise cette méthode exacte pour récupérer pour $D \in [[1, D_{m a x}]]$ la meilleure segmentation en $D$ segments minimisant le risque empirique i.e. :

\begin{matrix} {\hat{m}}_{D} = \underset{m \in ℳ (D)}{a r g m i n} ‖ Y - {\hat{μ}}_{m} ‖_{ℋ^{T}}^{2} . \end{matrix}

Nous présentons dans cette section l'algorithme de programmation dynamique appliquée au problème de détection de ruptures^[3]^,^[4]. Dans un premier temps, nous exprimons le risque empirique $‖ Y - {\hat{μ}}_{m} ‖_{ℋ^{T}}^{2}$ en fonction du noyau $k$ et de ${X_{i}}_{i \in [[1, T]]}$ à l'aide des deux résultats suivants ci-dessous.

On montre tout d'abord que, pour $m \in ℳ (D), d \in [[1, D - 1]]$ , pour $i \in [[τ_{d}, τ_{d + 1} - 1]]$ , ${({\hat{μ}}_{m})}_{i} = \frac{1}{τ_{d + 1} - τ_{d}} \sum_{j = τ_{d}}^{τ_{d + 1} - 1} Y_{j}$ .

Démonstration

Pour $u \in S_{m}^{ℋ}$ , notons, pour $λ \in m$ , $u_{λ}$ la valeur commune de $(u_{i})_{i \in λ}$ et, pour $g \in S_{m}^{ℋ}$ , ${\overline{g}}_{λ} = \frac{1}{c a r d (λ)} \sum_{i \in λ} g_{i}$ ,
Pour $g \in S_{m}^{ℋ}$ ,

\begin{matrix} ‖ u - g ‖_{ℋ^{T}}^{2} & = \sum_{λ \in m} \sum_{i \in λ} ‖ u_{λ} - g_{i} ‖_{ℋ}^{2} \\ = \sum_{λ \in m} \sum_{i \in λ} [‖ u_{λ} - {\overline{g}}_{λ} ‖_{ℋ}^{2} + ‖ g_{i} - {\overline{g}}_{λ} ‖_{ℋ}^{2} + 2 < u_{λ} - {\overline{g}}_{λ}, {\overline{g}}_{λ} - g_{i} >_{ℋ}] \\ = \sum_{λ \in m} [c a r d (λ) ‖ u_{λ} - {\overline{g}}_{λ} ‖_{ℋ}^{2}] + \sum_{λ \in m} \sum_{i \in λ} ‖ g_{i} - {\overline{g}}_{λ} ‖_{ℋ}^{2} + 2 \sum_{λ \in m} < u_{λ} - {\overline{g}}_{λ}, \sum_{i \in λ} ({\overline{g}}_{λ} - g_{i}) >_{ℋ} \\ = \sum_{λ \in m} [c a r d (λ) ‖ u_{λ} - {\overline{g}}_{λ} ‖_{ℋ}^{2}] + \sum_{λ \in m} \sum_{i \in λ} ‖ g_{i} - {\overline{g}}_{λ} ‖_{ℋ}^{2}, \end{matrix}

Ainsi, $‖ u - g ‖_{ℋ^{T}}^{2}$ est minimale si et seulement si $u_{λ} = {\overline{g}}_{λ}$ pour chaque $λ \in m$ .

On montre que, pour $m \in ℳ (D)$ ,

\begin{matrix} ‖ Y - {\hat{μ}}_{m} ‖_{ℋ^{T}}^{2} = \sum_{d \in [[1, D - 1]]} C_{τ_{d}, τ_{d + 1}}, \end{matrix}

avec, pour $(τ, τ^{'}) \in [[1, T]]^{2}, τ < τ^{'}$ , $C_{τ, τ^{'}} = \sum_{i = τ}^{τ^{'} - 1} k (X_{i}, X_{i}) - \frac{1}{τ^{'} - τ} \sum_{i = τ}^{τ^{'} - 1} \sum_{j = τ}^{τ^{'} - 1} k (X_{i}, X_{j}) .$

Démonstration

\begin{matrix} ‖ Y - {\hat{μ}}_{m} ‖_{ℋ^{T}}^{2} & = \sum_{d \in [[1, D - 1]]} \sum_{i \in [[τ_{d}, τ_{d + 1} - 1]]} ‖ Y_{i} - {({\hat{μ}}_{m})}_{i} ‖_{ℋ}^{2} \\ = \sum_{d \in [[1, D - 1]]} \sum_{i \in [[τ_{d}, τ_{d + 1} - 1]]} [‖ Y_{i} ‖_{ℋ}^{2} + ‖ \frac{1}{τ_{d + 1} - τ_{d}} \sum_{j = τ_{d}}^{τ_{d + 1} - 1} Y_{j} ‖_{ℋ}^{2} - 2 < Y_{i}, \frac{1}{τ_{d + 1} - τ_{d}} \sum_{j = τ_{d}}^{τ_{d + 1} - 1} Y_{j} >_{ℋ}] \\ = \sum_{d \in [[1, D - 1]]} (\sum_{i \in [[τ_{d}, τ_{d + 1} - 1]]} ‖ Y_{i} ‖_{ℋ}^{2} - \frac{1}{τ_{d + 1} - τ_{d}} \sum_{i = τ_{d}}^{τ_{d + 1} - 1} \sum_{j = τ_{d}}^{τ_{d + 1} - 1} < Y_{i}, Y_{j} >_{ℋ}) \\ = \sum_{d \in [[1, D - 1]]} (\sum_{i = τ_{d}}^{τ_{d + 1} - 1} k (X_{i}, X_{i}) - \frac{1}{τ_{d + 1} - τ_{d}} \sum_{i = τ_{d}}^{τ_{d + 1} - 1} \sum_{j = τ_{d}}^{τ_{d + 1} - 1} k (X_{i}, X_{j})) \\ = \sum_{d \in [[1, D - 1]]} C_{τ_{d}, τ_{d + 1}} . \end{matrix}

KPGD est l'implémentation du principe de la programmation dynamique à noyau appliquée au problème de détection de ruptures. Elle prend en paramètre la matrice de coût $C = (C_{i, j})_{(i, j) \in [[1, T]]} (i < j)$ et elle renvoie ${{\hat{m}}_{D}}_{D \in [[1, D_{m a x}]]}$ .

 algorithme KPGD ( $C$ ) 
   for  $d \in [[2, D_{m a x}]]$  do
     for  $τ^{'} \in [[d, T]]$  do
        $L_{d, τ^{'}} = min_{τ \leq τ^{'}} {L_{d - 1, τ} + C_{τ, τ^{'} + 1}}$ 
        $m_{d, τ^{'} + 1} = \underset{τ \leq τ^{'}}{a r g m i n} {L_{d - 1, τ} + C_{τ, τ^{'} + 1}}$ 
     end for        
   end for
 Fin algorithme

avec $\forall τ > 1, L_{1, τ} = C_{1, τ}$ .

La sélection de modèle^[4] permet de récupérer $\hat{D}$ un estimateur du nombre de segments $D^{*}$ . $\hat{D}$ est défini par

\begin{matrix} \hat{D} = \underset{D \in [[1, D_{m a x}]]}{a r g m i n} {‖ Y - {\hat{μ}}_{{\hat{m}}_{D}} ‖_{ℋ^{T}}^{2} + p e n ({\hat{m}}_{D})}, \end{matrix}

avec $p e n ({\hat{m}}_{D}) = \frac{C D}{n} (1 + \log (\frac{n}{D}))$ . La méthode utilisée pour calibrer la constante $C$ est l'heuristique de pente^[5]. On obtient ainsi un estimateur ${\hat{m}}_{\hat{D}}$ .

Segmentation binaire

L'heuristique de segmentation binaire^[6] est une méthode, fonctionnant par dichotomie, permettant de récupérer un minimiseur local ${\hat{m}}_{D}$ du risque empirique $‖ Y - {\hat{μ}}_{m} ‖_{ℋ^{T}}^{2}$ . Plus précisément, la segmentation binaire cherche à la première itération l'indice $j_{1} \in [[2, T - 1]]$ de l'instant de ruptures candidat qui minimise le risque empirique $‖ Y - {\hat{μ}}_{m} ‖_{ℋ^{T}}^{2}$ : cet indice est l'indice de notre premier instant de ruptures estimé. Puis, elle détermine récursivement, à la deuxième itération, deux instants de ruptures candidats $j_{2}, j_{3}$ sur chacun des segments délimités par les instants de ruptures estimés. Elle retient comme second instant de ruptures estimé celui (parmi ces deux instants de ruptures candidats) qui minimisent le risque empirique. Puis, on procède de la même manière pour les itérations suivantes. Nous illustrons sur un exemple le fonctionnement de l'algorithme utilisant le principe de la segmentation binaire :

Étape 1 : A l'itération $k = 1$ , on cherche $j_{1}$ qui minimise le risque empirique $‖ Y - {\hat{μ}}_{m_{1}} ‖_{ℋ^{T}}^{2}$ avec $m_{1} = {1, j, T}$ pour $j \in [[2, T - 1]]$ avec :

\begin{matrix} ‖ Y - {\hat{μ}}_{m_{1}} ‖_{ℋ^{T}}^{2} = C_{1, j} + C_{j, T} . \end{matrix}

$j_{1}$ est notre premier instant de ruptures estimé noté ${\hat{τ}}_{1}$ .

Étape 2 : A l'itération on $k = 2$ , on cherche $j_{2}, j_{3}$ minimisant le risque empirique sur $[[2, {\hat{τ}}_{1} - 1]]$ et $[[{\hat{τ}}_{1} + 1, T - 1]]$ respectivement. Par exemple,

j_{2} = \underset{j \in [[2, {\hat{τ}}_{1} - 1]]}{a r g m i n} ‖ Y_{1}^{{\hat{τ}}_{1}} - {\hat{μ}}_{m_{2}} ‖_{ℋ^{{\hat{τ}}_{1}}}^{2},

avec $m_{2} = (1, j, {\hat{τ}}_{1}, T)$ . Puis on choisit parmi les instants de ruptures candidats $j_{2}, j_{3}$ celui qui minimise le risque empirique $‖ Y - {\hat{μ}}_{m} ‖_{ℋ^{T}}^{2}$
et on le note ${\hat{τ}}_{2}$ . Puis, on continue récursivement.

Ainsi, au bout de $k$ itérations, on récupère une segmentation en $D = k + 1$ segments vérifiant :

\begin{matrix} {\hat{m}}_{D} = \underset{m \in ℳ_{D} ({\hat{τ}}_{1}, \dots, {\hat{τ}}_{D - 2})}{a r g m i n} ‖ Y - {\hat{μ}}_{m} ‖_{ℋ^{T}}^{2}, \end{matrix}

avec $ℳ_{D} ({\hat{τ}}_{1}, \dots, {\hat{τ}}_{D - 2})$ l'espace des segmentations en $D$ segments où les instants de ruptures estimés ${\hat{τ}}_{1}, \dots, {\hat{τ}}_{D - 2}$ ont été calculés aux itérations précédentes.

Une méthode alternative de segmentation binaire avec temps d'arrêt^[7] permet d'estimer le nombre de segments et donc de récupérer un estimateur ${\hat{m}}_{\hat{D}}$ de $m^{*}$ .

Notes et références

Articles

↑ [PDF] Alexandre Lung-Yut-Fong, « Détection de ruptures pour les signaux multidimensionnels. Application à la détection d’anomalies dans les réseaux. », HAL,‎ 2011 (lire en ligne)
↑ [PDF] Lucien Birgé et Pascal Massart, « Gaussian model selection. », Journal of the European Mathematical Society,‎ 1999, p. 203-268
↑ [PDF] Zaid Harchaoui et Olivier Cappé, « Retrospective Multiple Change-Point Estimation with Kernels. », IEEE,‎ 2007, p. 768-772 (lire en ligne)
1 2 [PDF] Sylvain Arlot, Alain Celisse et Zaid Harchaoui, « A kernel multiple change-point algorithm via model selection. », Arxiv,‎ 2016 (lire en ligne)
↑ Lucien Birgé et Pascal Massart, « Minimal Penalties for Gaussian Model Selection », Probability Theory and Related Fields, vol. 138,‎ 2007, p. 33--73 (lire en ligne).
↑ [PDF] Rebecca Killick, « Optimal detection of changepoints with a linear computational cost. », Journal of the American Statistical Association,‎ 2012, p. 1590--1598
↑ [PDF] Piotr Fryzlewicz, « Wild binary segmentation for multiple change-point detection. », The Annals of Statistics, vol. 42, n^o 6,‎ 2014, p. 2243--2281 (lire en ligne)

Portail des probabilités et de la statistique

[1] [PDF] Alexandre Lung-Yut-Fong, « Détection de ruptures pour les signaux multidimensionnels. Application à la détection d’anomalies dans les réseaux. », HAL,‎ 2011 (lire en ligne)

[2] [PDF] Lucien Birgé et Pascal Massart, « Gaussian model selection. », Journal of the European Mathematical Society,‎ 1999, p. 203-268

[3] [PDF] Zaid Harchaoui et Olivier Cappé, « Retrospective Multiple Change-Point Estimation with Kernels. », IEEE,‎ 2007, p. 768-772 (lire en ligne)

[kernel_multime_changepoint-4] 1 2 [PDF] Sylvain Arlot, Alain Celisse et Zaid Harchaoui, « A kernel multiple change-point algorithm via model selection. », Arxiv,‎ 2016 (lire en ligne)

[5] Lucien Birgé et Pascal Massart, « Minimal Penalties for Gaussian Model Selection », Probability Theory and Related Fields, vol. 138,‎ 2007, p. 33--73 (lire en ligne).

[6] [PDF] Rebecca Killick, « Optimal detection of changepoints with a linear computational cost. », Journal of the American Statistical Association,‎ 2012, p. 1590--1598

[7] [PDF] Piotr Fryzlewicz, « Wild binary segmentation for multiple change-point detection. », The Annals of Statistics, vol. 42, n^o 6,‎ 2014, p. 2243--2281 (lire en ligne)

[1]

[2]

[3]

[4]

[5]

[6]

[7]