\documentclass[XUPS,XML,SOM,Unicode,francais, NoFloatCountersInSection]{cedram}
\setcounter{tocdepth}{2}

\usepackage{../xups24}
\graphicspath{{xups24-01_figures}}

\let\preceq\preccurlyeq
\newcommand\mto{\mathchoice{\longmapsto}{\mapsto}{\mapsto}{\mapsto}}

%\XUPScorrections

\begin{document}
\frontmatter

\title[Introduction à la théorie de la persistance]{Introduction à la théorie de\nobreakspace la\nobreakspace persistance à travers un\nobreakspace exemple d'application}

\author{\firstname{Steve} \lastname{Oudot}}
\address{Inria Saclay - Ile-de-France \& École polytechnique,
Bât. Alan Turing,
1 rue Honoré d'Estienne d'Orves,
91120 Palaiseau}
\email{steve.oudot@inria.fr}
\urladdr{https://geometrica.saclay.inria.fr/team/Steve.Oudot/}

%\thanks{Journées X-UPS 2024. Introduction à l'analyse de données. Éditions de l'École polytechnique, 2024}

\begin{abstract}
Dans ce texte nous présentons de manière informelle les idées qui sous-tendent la théorie de la persistance, en nous appuyant sur un exemple particulier d’application de cette théorie au regroupement de données.\end{abstract}
\maketitle

\vspace*{-\baselineskip}
\tableofcontents
\mainmatter

Le but de ce texte est d'introduire quelques unes des principales idées sur lesquelles repose la théorie de la persistance, qui est le fondement mathématique de l'analyse topologique de données et le sujet central de cet ouvrage. Pour ce faire, nous allons nous placer dans un cadre restreint, celui de la persistance des pics d'une fonction réelle, et nous allons prendre pour prétexte l'application de la théorie à la classification non supervisée (aussi appelée \emph{regroupement}) en apprentissage machine. L'exposé alternera donc entre des sections de présentation du contexte applicatif et des sections plus formelles intro\-duisant les notions mathématiques. Parmi ces sections, la seule qui soit vraiment utile pour les textes suivants est la section~\ref{sec:pers_deg_0}, qui introduit la persistance dans notre cadre. Le lecteur qui ne s'intéresserait pas au contexte applicatif peut sans risque limiter sa lecture à cette seule section avant de passer au texte suivant.

\begin{figure}[htb]
\centering
\includegraphics[width=1\textwidth]{clustering}
\caption{Un échantillon de points dans~$\R^2$ (à gauche) et un regroupement de ces points (à droite).}
\label{fig:clustering}
\end{figure}

\section{Le problème du regroupement}
\label{sec:clustering}

Soit $P=\{p_1,\dots, p_n\}$ un ensemble fini de points de l'espace euclidien $\R^d$ -- dans le jargon on parle d'\emph{échantillon de points}, de \emph{jeu de données} ou encore de \emph{nuage de points}, selon le domaine scientifique considéré (statistiques, apprentissage machine, ou géométrie). On suppose que les points proviennent de $m$ classes distinctes, que l'on ne connaît pas -- $m$ lui-même n'est généralement pas connu non plus. L'objectif est donc (au besoin) de déterminer le nombre~$m$ de classes, puis de partitionner le nuage~$P$ en $m$ sous-ensembles appelés \emph{clusters}, $P=\bigsqcup_{\ell=1}^m C_\ell$, de manière à respecter au mieux les $m$ classes sous-jacentes. Voir la figure~\ref{fig:clustering} pour une illustration. Cette partition du nuage~$P$ s'appelle un \emph{regroupement} des points (\emph{clustering} en anglais).

Tel que formulé, le problème est clairement mal posé puisque, en l'absence d'informations complémentaires sur la manière dont les points sont obtenus à partir des classes sous-jacentes, il n'y a \hbox{aucun} moyen de déterminer si un regroupement du nuage est meilleur qu'un autre. Il faut donc formuler des \emph{hypothèses sur la génération des données}.

\begin{figure}[htb]
\centering
\includegraphics[width=1\textwidth]{mode-seeking}
\caption{À gauche: le nuage de points~$P$ de la figure~\ref{fig:clustering} superposé à la densité de probabilité~$f$ selon laquelle les points de~$P$ ont été échantillonnés iid. Les valeurs de~$f$ sont représentées par une palette de couleurs, du bleu (densité faible) à orange (densité forte). À droite: les six pics de~$f$ et leurs variétés stables (de la même couleur que le pic correspondant), ainsi qu'une représentation synthétique du flot de gradient de~$f$.}
\label{fig:mode-seeking}
\end{figure}

Chaque méthode de clustering vient avec son propre jeu d'hypothèses, duquel découle un critère de qualité sur les partitions de~$P$. Ici,~on va s'intéresser aux approches dites \emph{basées sur la densité}, qui formulent l'hypothèse suivante qu'on adoptera tout au long de ce texte -- voir la figure~\ref{fig:mode-seeking}~(gauche) pour une illustration:
\begin{hypothese}\label{hyp:sampling}
Les points $p_1,\dots, p_n$ sont échantillonnés de manière iid selon une mesure de probabilité~$\mu$, de densité~$f$ par rapport à la mesure de Lebesgue sur~$\R^d$. La mesure~$\mu$ comme sa densité~$f$ sont supposées \emph{inconnues}.
\end{hypothese}
Les classes sont alors définies comme des sous-ensembles deux à deux disjoints de~$\R^d$ associés aux maxima locaux -- aussi appelés \emph{modes} ou \emph{pics} -- de la densité~$f$. Plus précisément, elles correspondent aux \emph{variétés stables} des modes, c'est-à-dire qu'il y a exactement une classe par mode~$x$, et cette classe est définie comme étant le lieu des points de~$\R^d$ qui convergent asymptotiquement vers~$x$ lorsqu'ils sont poussés continûment le long du flot de gradient de la fonction~$f$ -- voir la figure~\ref{fig:mode-seeking}~(droite) pour une illustration. Nous allons maintenant formaliser cette notion en utilisant un peu de théorie de Morse, dont une référence classique est le livre de Milnor~\cite{milnor1963morse}.

\section{Définition mathématique des classes à partir de\nobreakspace la\nobreakspace densité\nobreakspace \texorpdfstring{$f$}{f}}
\label{sec:clustering_classes}

Comme on vient de le dire, pour définir les classes associées aux modes de notre densité~$f$ nous allons pousser les points de l'espace~$\R^d$ continûment le long du flot de gradient de~$f$. De toute évidence il nous faut faire des hypothèses de régularité sur~$f$ afin d'avoir un gradient et un flot bien définis, ainsi qu'un contrôle sur la convergence des trajectoires des points. Voici un jeu simplifié d'hypothèses, qu'on supposera vérifiées dans toute la suite du texte:
\begin{hypothese}\label{hyp:density}
On suppose que la densité~$f$:
\begin{enumeratei}
\item\label{item:i}
est lisse, de classe~$C^\infty$;
\item\label{item:ii}
s'annule à l'infini, c'est-à-dire que $\lim_{\|x\|\to\infty} f(x) = 0$;
\item\label{item:iii}
est \emph{de type Morse}, c'est-à-dire que ses \emph{points critiques}, \ie les points $x$ où le gradient $\nabla f(x)$ s'annule, sont en nombre fini et \emph{non dégénérés}, c'est-à-dire que la matrice hessienne de $f$ en chacun de ces points est inversible.
\end{enumeratei}
\end{hypothese}
L'hypothèse~\eqref{item:iii} joue un rôle clé dans la suite. Elle implique en particulier que les points critiques de~$f$ sont isolés dans~$\R^d$. Bien qu'elle puisse paraître plus restrictive que les autres hypothèses a priori, puisqu'elle impose des conditions sur les quantités différentielles d'ordre~1 et~2 de~$f$, il s'avère que les fonctions de type Morse forment un ouvert dense pour la topologie~$C^2$ dans l'espace des fonctions~$C^\infty$, donc on ne perd presque rien en généralité en ajoutant l'hypothèse~\eqref{item:iii}.

Comme la fonction~$f$ est de classe~$C^\infty$, son champ de gradient est localement lipschitzien et peut donc être intégré en un flot continu $\Phi\colon \R^+\times \R^d\to\R^d$. Plus précisément, par le théorème de Cauchy-Lipschitz global, la ligne de flot $\varphi_x\colon t\!\mto\! \Phi(t,x)$ issue d'un point~\hbox{$x\!\in\!\R^d$} est l'unique solution de l'équation différentielle ordinaire suivante:
\[
\begin{cases}
\dot{\varphi}_x(t) = \nabla f(\varphi_x(t))\\
\varphi_x(0) = x
\end{cases}
\]
Cette solution dépend continûment à la fois du paramètre~$t$ et de la condition initiale~$x$, d'où la continuité du flot~$\Phi$.

On regarde maintenant, pour tout point~$x\in\R^d$, si et où converge la ligne de flot~$\varphi_x$ lorsque $t\to +\infty$. Le cas particulier où $f(x)=0$ est trivial: en tant que densité de probabilité, $f$ est positive ou nulle, donc $x$ est un minimum local de~$f$ et, à ce titre, lui-même un point critique de~$f$, stationnaire pour le flot~$\Phi$ par définition. Pour le cas général $f(x)>0$, l'hypothèse que $f$ s'annule à l'infini (sur $\R^d$ qui est localement compact) implique que l'ensemble de sur-niveau
\[
\{y\in\R^d \mid f(y)\geq f(x)>0\},
\]
dans lequel se trouve l'image de la ligne de flot~$\varphi_x$, est compact; il~s'ensuit alors que la ligne de flot $\varphi_x$ converge bien lorsque $t\to+\infty$, et par définition la limite est un point critique de~$f$. Ainsi, à la limite, le flot amène tous les points de~$\R^d$ à des points critiques de~$f$.

\begin{definition}
La \emph{variété stable} d'un point critique~$x$ de~$f$ est l'ensemble des points $y\in\R^d$ tels que $x=\lim_{t\to +\infty}\varphi_y(t)$.
\end{definition}

Il découle de ce qui vient d'être dit que les variétés stables des points critiques partitionnent l'espace~$\R^d$. Toutefois, tous les points critiques ne sont pas des pics: il y a également les minima locaux, ainsi que les points critiques de type \emph{selle}. Dans la suite on ne retiendra que les variétés stables des pics.

\begin{definition}
Les classes correspondant à la densité~$f$ sont par définition les variétés stables des pics de~$f$.
\end{definition}

L'idée derrière le fait de ne regarder que les variétés stables des pics est que, en général, identifier les points critiques requiert d'évaluer des quantités différentielles, typiquement le gradient de~$f$ (voire la matrice hessienne si on veut déterminer le type des points critiques), dont le calcul est instable numériquement. Les pics, quant à eux, ne nécessitent pas de quantités différentielles pour être caractérisés (définition~\ref{def:local_max}), leur calcul en pratique s'avère donc beaucoup plus stable numériquement.

En contrepartie, on peut s'interroger sur la pertinence de laisser de côté les autres types de points critiques, impliquant de fait que l'union des classes ne couvre pas $\R^d$ tout entier. En réalité, la théorie de Morse nous garantit que les classes couvrent presque tout l'espace:

\begin{proposition}
Sous l'hypothèse~\ref{hyp:density}, le complémentaire de l'union des variétés stables des pics de~$f$ est une union finie de sous-variétés différentielles de~$\R^d$ de codimensions strictement positives, et donc de mesure de Lebesgue nulle dans~$\R^d$.
\end{proposition}

En pratique, la probabilité qu'un point~$p_i$ de notre échantillon~$P$ ne soit pas parmi les classes est donc nulle, puisque les $p_i$ sont échantillonnés selon la mesure $\mu$ qui a une densité par rapport à la mesure de Lebesgue. Ainsi, nos hypothèses initiales et notre définition des classes sont génériquement compatibles avec notre problème.

\section{Calcul des clusters à partir du nuage de points~\texorpdfstring{$P$}{P}}
\label{sec:clustering_computation}

Comme on l'a supposé dans l'hypothèse de départ~\ref{hyp:sampling}, la mesure~$\mu$ et sa densité~$f$ ne sont pas connues en pratique. Il nous va donc falloir trouver un moyen de simuler la montée de gradient à partir des données $p_1,\dots, p_n$ pour pouvoir former des clusters qui approximent les classes. Il existe tout un éventail de méthodes pour ce faire. Les unes, comme par exemple \emph{mean-shift}~\cite{cm-ms-02}, adoptent une approche purement numérique en tentant d'approximer localement le gradient de~$f$ puis de simuler la montée de gradient continue par une montée de gradient approximé discrète dans~$\R^d$. Les autres, plus anciennes comme celle que nous allons voir ici~\cite{knf-gtanp-76}, adoptent une approche combinatoire en remplaçant $\R^d$ par un objet discret appelé \emph{graphe de voisinage}, construit à partir des données, dans lequel le gradient est approximé en chaque sommet par une arête incidente.

\subsubsection*{Pré-traitement}

En pré-traitement on construit le graphe de voisinage et on approxime la densité aux sommets du graphe.

\begin{figure}[htb]
\centering
\includegraphics[width=1\textwidth]{Koontz}
\caption{À gauche: le graphe~$G$ de $r$-voisinage sur le nuage de points~$P$ de la figure~\ref{fig:clustering}, pour la valeur de~$r$ montrée au centre. Les arêtes du graphe sont en noir, tandis que les valeurs de l'estimateur~$\hat f$ aux sommets sont représentées comme à la figure~\ref{fig:mode-seeking}. À droite: le résultat de l'algorithme sur cette entrée, avec une couleur distincte par cluster. Les arêtes représentées sur le dessin sont celles des arbres de la forêt couvrante de~$G$ calculée par l'algorithme.}
\label{fig:Koontz}
\end{figure}

\begin{definition}
Un \emph{prédicat de voisinage} est une fonction symétrique $P\times P\to \{0,1\}$ qui vaut~$0$ sur la diagonale $\{(p_i,p_i) \mid p_i\in P\}$.
\end{definition}

\begin{definition}
Étant donné un prédicat de voisinage~$\sigma\colon P\times P\to \{0,1\}$, le \emph{graphe de voisinage} correspondant est le graphe combinatoire $G=(P,E)$ non orienté dont les sommets sont les points de~$P$ et les arêtes forment l'ensemble $E=\{(p_i, p_j)\in P\times P \mid \sigma(p_i, p_j)=1\}$.
\end{definition}

\begin{exemple}
Voici deux constructions classiques de graphes de voisinage, dont la première est illustrée dans la figure~\ref{fig:Koontz}~(gauche):
\begin{itemize}
\item le graphe de \emph{$r$-voisinage}, pour un réel~$r\geq 0$, correspond au prédicat $(p_i, p_j) \mto \one_{p_i\neq p_j}\, \one_{\|p_i-p_j\|_2\leq r}$ qui teste si deux points donnés sont à distance euclidienne au plus~$r$ l'un de l'autre;
\item le graphe de \emph{$k$-voisinage}, pour un entier $k\geq 1$, correspond au prédicat $(p_i, p_j) \mto \one_{p_i\neq p_j}\, \one_{p_j\in \text{ppv}_k(p_i)}\, \one_{p_i\in \text{ppv}_k(p_j)}$, où $\text{ppv}_k(p_\ell)$ dési\-gne l'ensemble des $k$ plus proches voisins de~$p_\ell$ parmi les points du nuage~$P$ pour la distance euclidienne.
\end{itemize}
\end{exemple}

Une fois un tel graphe de voisinage~$G=(P,E)$ construit à partir de~$P$, on calcule une estimation $\hat f(p_i)$ de la densité~$f$ en chaque point $p_i\in P$. L'estimation de la densité est en soi un sujet à part entière, qui sort du cadre de cet exposé. Notons simplement que le domaine des statistiques nous fournit tout un éventail d'estimateurs ayant chacun des propriétés spécifiques. Ici nous ferons simplement l'hypothèse (relativement forte mais classique) que l'estimateur~$\hat f$ appro\-xime la densité~$f$ en norme sup sur~$P$, c'est-à-dire que l'on peut borner l'erreur de l'estimateur comme suit:
\begin{equation}\label{eq:density_error}
\|\hat f - f\|_\infty \defeq \max_{1\leq i\leq n} |\hat f(p_i)-f(p_i)| \leq \e(n)
\end{equation}
où la quantité~$\e(n)$ dépend uniquement de $n$, pas des points du nuage~$P$, et tend vers~$0$ lorsque $n\to +\infty$. La borne $\e(n)$ en elle-même n'a pas d'importance ici, et pour simplifier nous la supposerons vérifiée de manière déterministe et non pas seulement avec forte probabilité comme c'est le cas en pratique. Ainsi nous avons équipé notre graphe de voisinage $G=(P,E)$ d'un champ scalaire $\hat f\colon P\to\R^+$ appro\-ximant la densité~$f$ aux sommets -- voir encore la figure~\ref{fig:Koontz}~(gauche). C'est l'entrée que nous fournissons à l'algorithme.

\subsubsection*{L'algorithme} En chaque sommet $p_i$ du graphe nous approximons le gradient de~$f$ par l'arête reliant $p_i$ à son voisin $p_j$ dont l'estimée~$\hat f(p_j)$ est la plus élevée, \emph{à condition que celle-ci soit plus élevée que~$\hat f(p_i)$}. Dans le cas contraire, $p_i$ est un maximum local de~$\hat f$ dans le graphe~$G$ et on le déclare donc comme étant un pic de la densité.

L'ensemble des arêtes ainsi sélectionnées forme une \emph{forêt couvrante} de~$G$, c'est-à-dire un ensemble de sous-graphes qui sont des arbres (\ie des sous-graphes connexes sans cycles) et qui couvrent des sous-ensembles deux à deux disjoints de sommets dont l'union est le nuage~$P$ tout entier\footnote{Techniquement, certains sommets de~$G$ peuvent être isolés, auquel cas on les ajoute à la forêt en tant qu'arbres singletons.}. Ces arbres sont les clusters produits par l'algorithme. Chacun contient une unique racine, son sommet dont la valeur de $\hat f$ est la plus élevée, qui par construction est un maximum local de $\hat f$ dans~$G$ et sert donc d'approximation pour un éventuel pic de~$f$. L'arbre en lui-même sert d'approximation pour la variété stable associée à ce pic dans~$\R^d$.

\begin{figure}[htb]
\centering
\includegraphics[width=0.5\textwidth]{density_noise}
\caption{Graphe de l'estimateur de densité~$\hat f$ restreint au nuage de points~$P$ de la figure~\ref{fig:clustering}.}
\label{fig:density_noise}
\end{figure}

\subsubsection*{Résultat} La figure~\ref{fig:Koontz}~(droite) montre un exemple de résultat de l'algorithme, qui, il faut l'avouer, est de piètre qualité. Ce qui frappe immé\-diatement, c'est la multiplication des clusters (plusieurs dizaines) par rapport au nombre de classes sous-jacentes (six seulement). Ceci s'explique essentiellement par le fait que l'estimateur~$\hat f$ est beaucoup plus bruité que la densité~$f$, n'étant qu'une approximation en norme sup d'après~\eqref{eq:density_error} -- voir la figure~\ref{fig:density_noise} pour une illustration du bruit dans l'estimateur. Ainsi, en plus de quelques pics \og légitimes\fg associés aux pics de~$f$ dans~$\R^d$, s'ajoute dans le graphe~$G$ tout un tas de pics fallacieux dus au bruit dans~$\hat f$. Afin de distinguer parmi les pics ceux qui sont légitimes de ceux qui ne le sont pas, nous allons utiliser une notion mathématique qui quantifie l'importance de chaque pic: la~\emph{persistance}. Cette notion va également nous fournir une \emph{hiérarchie sur les pics} de~$\hat f$ dans~$G$, hiérarchie qui va nous permettre de réparer le clustering produit par l'algorithme en fusionnant les clusters associés aux pics fallacieux avec les clusters de leurs parents dans la hiérarchie.

\section{La persistance des pics d'une fonction réelle}
\label{sec:pers_deg_0}

Alors que le terme \emph{persistance} est utilisé communément en analyse topologique de données, dans le cas particulier des pics il renvoie à un concept plus ancien, la \emph{proéminence}, introduit par les alpinistes et les géographes, que nous allons présenter en premier afin de nourrir l'intuition du lecteur.

\subsection{Définition des alpinistes et géographes}
\label{sec:def_geographes}

\begin{figure}[htb]
\centering
\includegraphics[width=0.7\textwidth]{proeminence_geo}
\caption{Les flèches verticales montrent la proéminence de trois pics sur une île. Une ligne pointillée horizontale relie chaque pic (excepté le plus haut) à son col le plus élevé.\\{\footnotesize Source: \emph{Wikipedia} (\url{https://en.wikipedia.org/wiki/File:Relative-height.png}, image sous licence CC BY-SA 3.0 Deed).}}
\label{fig:proeminence_geo}
\end{figure}

Avec\footnote{Le contenu de cette sous-section est largement repris de l'article correspondant sur \emph{Wikipedia} (\href{https://fr.wikipedia.org/wiki/Pro\%C3\%A9minence}{\texttt{fr.wikipedia.org/wiki/Proéminence}}).}
 le développement de l'exploration alpine dans la deuxième moitié du \textsc{xix}\ieme siècle, des listes de sommets, conquis ou à conquérir, ont commencé à émerger. Rapidement s'est posée la question de déterminer ce qu'est un sommet. En effet, sans critère restrictif, n'importe quelle antécime, épaule, ou même pierre pourrait être vue comme un sommet en elle-même. Le critère principal retenu pour distinguer les sommets des autres types de protubérances a été celui de la proéminence, qui doit être supérieure à un certain seuil pour que la protubérance puisse être considérée comme un sommet indépendant. Ce seuil varie d'un classement à l'autre: de 30~mètres (la longueur de corde de l'alpinisme classique) pour la liste officielle des 82~sommets des Alpes de plus de 4000~mètres, à 1500~mètres pour les sommets dits \emph{ultra-proéminents}.

La proéminence pour les alpinistes et les géographes se définit comme \og la différence d'altitude entre un sommet donné et l'ensellement ou le col le plus élevé permettant d'atteindre une cime encore plus haute\fg. Autrement dit, c'est \og le dénivelé minimum de la descente à parcourir pour pouvoir remonter sur un sommet plus élevé\fg. Voir la figure~\ref{fig:proeminence_geo} pour une illustration. Notons que l'ensellement peut se situer au niveau de la mer mais pas au-dessous. Ainsi, le plus haut pic d'une île a une proéminence égale à sa hauteur. Les tables~\ref{tab:altitude} et~\ref{tab:proeminence} donnent les listes des dix sommets les plus hauts du monde d'une part, des dix sommets les plus proéminents du monde d'autre part, et comme on peut le constater elles sont bien différentes.

\begin{table}[t]
\caption{Liste des 10 sommets les plus hauts du monde, classés par ordre décroissant d'altitude. Le seuil minimal de proéminence requis ici est de 500~mètres.}
\label{tab:altitude}
\centering
\smaller
\begin{tabular}{|l|l|r|r|}
\hline
sommet & continent & altitude~(m) & proéminence~(m) \\
\hline
\hline
Everest & Asie & 8849 & 8849 \\
K2 & Asie & 8611 & 4020 \\
Kangchenjunga & Asie & 8586 & 3922 \\
Lhotse & Asie & 8516 & 610 \\
Makalu & Asie & 8485 & 2378 \\
Cho Oyu & Asie & 8188 & 2340 \\
Dhaulagiri I & Asie & 8167 & 3357 \\
Manaslu & Asie & 8163 & 3092 \\
Nanga Parbat & Asie & 8126 & 4608 \\
Annapurna I & Asie & 8091 & 2984\\
\hline
\end{tabular}
\end{table}

\begin{table}[t]

\caption{Liste des 10 sommets les plus proéminents du monde, classés par ordre décroissant de proéminence.}
\smaller
\label{tab:proeminence}
\centering
\begin{tabular}{|l|l|r|r|}
\hline
sommet & continent & altitude~(m) & proéminence~(m) \\
\hline
\hline
Everest & Asie & 8849 & 8849 \\
Aconcagua & Amérique & 6960 & 6960 \\
Denali & Amérique & 6191 & 6155 \\
Kilimanjaro & Afrique & 5895 & 5885 \\
Pico Cristóbal Colón & Amérique & 5700 & 5509 \\
Mont Logan & Amérique & 5959 & 5250 \\
Pico de Orizaba & Amérique & 5636 & 4922 \\
Massif Vinson & Antarctique & 4892 & 4892 \\
Puncak Jaya & Océanie & 4884 & 4884 \\
Mont Elbrouz & Europe & 5642 & 4741 \\
\hline
\end{tabular}
\end{table}

\subsection{Définition mathématique}
\label{sec:def_mathematiciens}

On va maintenant définir formellement la proéminence, ou persistance. Pour cela on fixe un espace topologique~$X$ et une fonction $f\colon X\to\R$, sans plus d'hypothèses pour le moment -- des hypothèses sur~$f$ et~$X$ apparaîtront au cours de l'exposé. Notons dès à présent que l'approche adoptée ici pour définir la persistance transcrit directement les idées de la section~\ref{sec:def_geographes} en requérant peu d'outils mathématiques. En contrepartie, elle donne lieu à des énoncés parfois inutilement techniques et se généralise mal -- voir à ce propos l'hypothèse~\ref{hyp:CCs} et les exemples et commentaires qui l'entourent. La bonne approche, fondée sur l'homologie, sera adoptée dans le texte~\cite{chap-pers-1} -- voir en particulier la remarque~5.11.

\begin{definition}\label{def:local_max}
Un point $x\in X$ est un \emph{maximum local} (ou \emph{pic}) de $f$ s'il existe un voisinage $U$ de $x$ dans~$X$ tel que $f(x) = \max_{U} f$.
\end{definition}

Pour quantifier la proéminence des pics, on regarde l'évolution des composantes connexes par arc dans les \emph{sur-niveaux} de la fonction~$f$ alors que le niveau diminue progressivement de $+\infty$ jusqu'à $-\infty$.

\begin{definition}\label{def:sous-sur-niveau}
Étant donné un niveau $t\in\R$, l'\emph{(ensemble de) sur-niveau} de~$f$ associé est $f^{-1}([t, +\infty))$.
\end{definition}

\begin{definition}\label{def:max_cluster}
Étant donné un pic $x\in X$ de $f$, pour tout niveau $t\leq f(x)$ on définit $C_t(x)$ comme étant la composante connexe par arc du sur-niveau $f^{-1}([t,+\infty))$ à laquelle appartient~$x$.
\end{definition}

Lorsque $t$ diminue, le sur-niveau de~$f$ associé ne fait que croître, ainsi que ses composantes connexes par arc. En conséquence:
\begin{corollaire}\label{cor:Ct_to_Ct'}
Étant donné un pic $x\in X$ de $f$, pour tous niveaux $t\leq t'\leq f(x)$ on a $C_{t}(x)\supseteq C_{t'}(x)$.
\end{corollaire}
Ainsi, informellement, on peut traquer la croissance de la composante connexe par arc contenant notre pic~$x$ tandis que l'on abaisse le niveau~$t$, et détecter la première valeur de~$t$ à laquelle cette composante fusionne avec celle d'un pic plus élevé: à cette valeur $t=h(x)$ particulière, $x$ émerge comme un pic secondaire d'une montagne plus élevée, et sa persistance est donnée par la différence $f(x)-h(x)\geq 0$. Dans le cas particulier où la composante de $x$ fusionne avec celle d'un pic de même hauteur, on départage les deux pics en désignant l'un comme étant secondaire de l'autre de manière arbitraire.

Formellement, on suppose donné un ordre total $\preceq$ sur $X$ (ce qui en toute généralité requiert une version faible de l'axiome du choix) et
on considère l'ordre lexicographique suivant sur~$X$:
\begin{equation}\label{eq:pics_lexico}
y\geq x \quad \Longleftrightarrow \quad \begin{cases} f(y)>f(x)& \text{ou} \\ f(y)=f(x) & \text{et}\quad y\preceq x \end{cases}
\end{equation}
On note~$>$ l'ordre total strict associé:
\[ y > x \quad \Longleftrightarrow \quad y\geq x\quad \text{et}\quad y\neq x \]
\begin{definition}\label{def:proeminence}
Pour tout pic $x\in X$ de $f$ on définit:
\begin{itemize}
\item l'instant de \emph{naissance} de $x$ comme étant la valeur $f(x)\in\R$;
\item l'instant de \emph{décès} de $x$ par
$h(x)=\sup\,I(x) \in\R\cup\{-\infty\}$, où~$I(x) = \{t\leq f(x) \mid \exists y>x\ \text{pic de~$f$ tel que}\ C_t(y)=C_t(x)\}$;
\item l'\emph{intervalle de persistance} de~$x$ comme étant l'intervalle semi-ouvert $]h(x), f(x)]\subseteq \R$, que l'on maintient ouvert à gauche par convention car cette extrémité peut être à l'infini;
\item la \emph{persistance} (ou \emph{proéminence} en géographie) de~$x$ comme étant la différence $f(x)-h(x)\in \R^+ \cup \{+\infty\}$.
\end{itemize}
\end{definition}

\begin{remarque}
Pour les pics $x$ tels que $h(x)=-\infty$, la proéminence telle que définie ici est infinie, tandis que pour les géographes la proéminence de ces pics est égale à leur hauteur (voir la section~\ref{sec:def_geographes} et en particulier la figure~\ref{fig:proeminence_geo}).
\end{remarque}

\begin{definition}\label{def:barcode_deg0}
Le \emph{code-barres} de $f$ est le multi-ensemble des intervalles de persistance des pics de~$f$. La \emph{multiplicité} d'un intervalle est le nombre (potentiellement infini) de ses occurrences dans le multi-ensemble.
\begin{figure}[htb]
\centering
\includegraphics[width=0.6\textwidth]{sinxcos3x}
\caption{Graphe (en brun) et code-barres (en violet) de la fonction $x\mto - \sin(x) \, \cos(3x)$ sur l'intervalle $X=[-5, 2]$ muni de l'ordre usuel sur les réels. Par souci de clarté, les inter\-valles de persistance sont ancrés à l'aplomb des pics correspondants de la fonction. Les liens de parenté formant la hiérarchie des pics sont matérialisés par des flèches courbes vertes reliant chaque pic à son parent lorsqu'il existe.}
\label{fig:-sinxcos3x}
\end{figure}
\end{definition}

\begin{exemple}\label{ex:-sinxcos3x}
Considérons la fonction $x\mto - \sin(x) \, \cos(3x)$ sur l'intervalle $X=[-5, 2]\subset\R$ muni de l'ordre usuel sur les réels, dont le graphe est représenté à la figure~\ref{fig:-sinxcos3x}. Cette fonction possède cinq pics, aux abscisses $x=-5$, $x=\arctan\sqrt{2+\sqrt{\sfrac{11}{3}}} + k\,\pi$ et \hbox{$x=-\arctan\sqrt{2-\sqrt{\sfrac{11}{3}}} + k\,\pi$} pour $k\in \{-1,\,0\}$. Les intervalles de son code-barres sont, de gauche à droite:
\[ ]-a,\, c] \quad\quad ]-b,\, b] \quad\quad ]-\infty,\, a] \quad\quad ]-b,\, b] \quad\quad ]-a,\, a] \]
où
\begin{align*}
a &= \left(3\sqrt{11/3}+5\right)\left(2+\sqrt{11/3}\right)^{1/2}\left(3+\sqrt{11/3}\right)^{-2} \ \approx\ 0.88\\
b &= \left(3\sqrt{11/3}-5\right)\left(2-\sqrt{11/3}\right)^{1/2}\left(3-\sqrt{11/3}\right)^{-2} \ \approx\ 0.19\\
c &= \sin(5) \cos(15)\ \approx\ 0.73
\end{align*}
\end{exemple}

Les notions d'intervalle de persistance et de code-barres sont défi\-nies pour toute fonction~$f\colon X\to\R$, toutefois elles n'ont de sens que lorsqu'on fait l'hypothèse que les composantes connexes par arc contenant les pics couvrent l'intégralité des sur-niveaux de la fonction, c'est-à-dire:
\begin{hypothese}\label{hyp:CCs}
\(\displaystyle \forall t\in \R, \quad f^{-1}([t, +\infty[) =\bigcup_{\begin{smallmatrix}x\ \mathrm{pic\ de}\ f\\f(x)\geq t\end{smallmatrix}} C_t(x) \)
\end{hypothese}
Dans le cas contraire en effet, des composantes connexes par arc dans les sur-niveaux peuvent être ignorées et donc des barres ne pas apparaître ou bien apparaître avec des extrémités erronées dans le code-barres, comme dans les exemples ci-dessous:

\begin{exemple}\label{ex:pic-pers_bad}
La fonction identité sur $\R$ ou sur l'intervalle $]0,1[$ n'a aucun pic et donc son code-barres est vide. De même pour la fonction
\[
x\mto \begin{cases}
1-|x|&\text{si}\ x\neq 0\\0&\text{si}\ x=0
\end{cases}\quad \text{sur l'intervalle $[-1,1]$}.
\]
\end{exemple}

\begin{exemple}
La fonction $f\colon x\mto -(x^3+2) \exp(-x)$ sur $\R^+$ a un unique pic en $x=1$, dont la proéminence est $+\infty$ alors qu'on s'attendrait à ce qu'elle soit finie car $f(1) = -\sfrac{3}{e}$ et $\lim_{t\to +\infty} f= 0$.
\end{exemple}

À partir de maintenant nous supposerons donc l'hypothèse~\ref{hyp:CCs} véri\-fiée. C'est le cas par exemple lorsque $X$ est compact et $f$ est continue, ou encore lorsque $X=\R^d$ et $f$ est continue, positive ou nulle et s'annule à l'infini comme sous l'hypothèse~\ref{hyp:density}.

\subsubsection*{Hiérarchie des pics}
En plus du code-barres, nous pouvons définir une notion de parent et de là une hiérarchie sur les pics. Pour cela nous faisons l'hypothèse additionnelle suivante, également vérifiée sous l'hypo\-thèse~\ref{hyp:density}:
\begin{hypothese}\label{hyp:finite_peaks}
Le nombre de pics de la fonction~$f$ est fini.
\end{hypothese}

\begin{exercice}\label{lem:parent}
Soit $x\in X$ un pic de~$f$ dont la proéminence est finie ($h(x)>-\infty$). Montrer que l'ensemble $I(x)$ de la définition~\ref{def:proeminence} est alors non vide, de même que l'ensemble
\[ J(x)\ =\ \left\{ y\ \text{pic de~$f$} \mid y >x\ \text{et}\ C_t(y)=C_t(x)\ \forall t\in I(x)\right\} \]
et que ce dernier
admet un maximum pour l'ordre~$\geq$ de l'équation~\eqref{eq:pics_lexico}.
\end{exercice}

Le maximum de~$J(x)$ est appelé le \emph{parent} de~$x$. Il
est strictement supérieur à~$x$ pour l'ordre~$\geq$ donc la relation de parenté induit une hiérarchie sur les pics, au sommet de laquelle se trouvent les pics de proéminence infinie. La figure~\ref{fig:-sinxcos3x} montre cette hiérarchie pour la fonction de l'exemple~\ref{ex:-sinxcos3x}.

\begin{figure}[htb]
\centering
\includegraphics[width=0.75\textwidth]{barcode_diagram}
\caption{Le code-barres issu de la figure~\ref{fig:-sinxcos3x} (à gauche) et son diagramme de persistance correspondant (à droite). La~multiplicité du point $(b,-b)$ dans le diagramme est indiquée entre parenthèses.}
\label{fig:barcode_diagram}
\end{figure}

\begin{figure}[htb]
\centering
\includegraphics[width=1\textwidth]{bottleneck}
\caption{Deux fonctions $[0,1]\to\R$ (à gauche) et leurs diagrammes de persistance (à droite) avec, marqué par des segments en noir, l'appariement optimal entre les points des deux diagrammes pour la distance de transport. Dans cet exemple, la distance de transport est légèrement inférieure à la différence des deux fonctions en norme sup.}
\label{fig:bottleneck}
\end{figure}

\subsubsection*{Diagrammes de persistance et stabilité}
Une autre manière de représenter graphiquement les codes-barres est sous la forme de multi-ensembles de points dans le plan étendu $\R\times [-\infty, +\infty[$, appelés \emph{diagrammes de persistance}, dans lesquels chaque copie du point $(a,b)$ correspond à une copie de l'intervalle $]b,a]$ dans le code-barres asso\-cié, avec $a > b \geq -\infty$. Voir la figure~\ref{fig:barcode_diagram} pour une illustration. Cette représentation alternative contient la même information mais elle offre l'avantage de montrer les codes-barres comme des nuages de points, ou plutôt comme des mesures empiriques, plus facilement interprétables.
Par ailleurs, la métrique naturelle entre codes-barres, appelée \emph{distance du goulot de bouteille} (\emph{bottleneck distance} en \hbox{anglais}), peut être interprétée comme une distance de transport partiel entre diagrammes de persistance, eux-mêmes vus comme des mesures empiriques. La définition précise de cette distance sera donnée au texte~\cite{sec-stabilite}, mais en attendant nous pouvons déjà la visualiser sur l'exemple de la figure~\ref{fig:bottleneck}. L'interprétation des deux diagrammes dans l'exemple est que chacune des deux fonctions considérées possède trois pics proéminents, le reste des pics (dont les points correspondants dans les diagrammes sont localisés près de la diagonale $y=x$) pouvant être considéré comme du bruit. Le plan de transport optimal asso\-cié à la distance du goulot de bouteille entre les diagrammes donne un appariement explicite entre les pics proéminents des deux fonctions, et nous indique par ailleurs d'ignorer les pics non proéminents (en appariant avec la diagonale leurs points correspondants dans les diagrammes).

Ainsi, les diagrammes de persistance fournissent une représentation synthétique des intervalles de persistance des pics d'une fonction, tandis que la distance de transport entre diagrammes indique comment mettre en correspondance au mieux les pics de fonctions différentes selon leurs intervalles de persistance. Cette observation empirique est appuyée par un résultat fondamental de la théorie de la persistance, appelé le \emph{théorème de stabilité}, qui dit en substance que l'opérateur~$D$, qui associe à toute fonction $X\to\R$ son diagramme de persistance lorsque celui-ci existe, est 1-lipschitzien. Dans notre contexte le résul\-tat s'énonce de la manière suivante -- voir encore la figure~\ref{fig:bottleneck} pour une illustration:
\begin{theoreme}\label{thm:stability_deg0}
Pour toutes fonctions $f,g\colon X\to\R$ vérifiant les hypothèses~\ref{hyp:CCs} et~\ref{hyp:finite_peaks}, on a l'inégalité suivante, où $\bottleneckd$ désigne la distance du goulot de bouteille et $\|\cdot\|_\infty$ désigne la norme sup sur~$X$:
\[ \bottleneckd(D(f),\, D(g)) \leq \|f-g\|_\infty \]
\end{theoreme}
Dans le texte~\cite{sec-stabilite} nous énoncerons et démontrerons ce résultat dans le cadre général de la persistance des fonctions réelles, pas seulement celle de leurs pics.

\section{Retour à l'application au regroupement}
\label{sec:back_to_clustering}

Revenons maintenant à notre application et utilisons la persistance pour corriger les défauts de l'algorithme de clustering présenté à la section~\ref{sec:clustering_computation}. Nous n'allons en fait pas modifier l'algorithme en lui-même, mais plutôt ajouter un post-traitement des clusters qu'il produit.

\subsubsection*{Post-traitement}
Étant donné le regroupement $P=\bigsqcup_{\ell=1}^m C_\ell$ produit par l'algorithme sur le graphe de voisinage~$G=(P,E)$, ainsi que les pics $x_1,\dots, x_m\in P$ de~$\hat f$ dans~$G$ associés à chacun des clusters $C_1,\dots, C_m$, on calcule l'intervalle de persistance de chaque pic~$x_\ell$ à l'intérieur du graphe~$G$, ainsi que son parent (s'il existe) parmi les autres pics. Les détails du calcul importent peu, il suffit de savoir qu'il est possible de le faire en un temps quasi-linéaire en la taille du graphe, par un algorithme similaire à celui de Kruskal pour le calcul d'arbre couvrant minimal~\cite[section~23.2] {clrs-ia-09}.
On obtient ainsi le code-barres de $\hat f$, vue comme une fonction réelle sur le graphe $G$, vu lui-même comme un espace topologique stratifié par ses sommets et ses arêtes, avec interpolation linéaire des valeurs de $\hat f$ le long des arêtes. Voir la figure~\ref{fig:ToMATo_PD} pour une illustration.

\begin{figure}[htb]
\centering
\includegraphics[width=.9\textwidth]{ToMATo_PD}
\caption{À gauche, en rouge et noir: les pics de l'estimateur de densité~$\hat f$ dans le graphe de voisinage~$G$ de la figure~\ref{fig:Koontz}. À droite: le diagramme de persistance de~$\hat f$, avec, en rouge, six points qui se démarquent clairement du reste du diagramme. Ces points correspondent aux intervalles de persistance des six pics de~$\hat f$ en rouge sur la figure de gauche, que l'on peut mettre en correspondance avec les six pics de la vraie densité~$f$ (qui sont tout proches, voir la figure~\ref{fig:mode-seeking}), le reste des pics de~$\hat f$ (en noir) se répartissant dans des zones de faible norme du gradient de~$f$ et pouvant être considéré comme du bruit d'après le diagramme de persistance de~$\hat f$.}
\label{fig:ToMATo_PD}
\vspace*{-\baselineskip}\enlargethispage{.5\baselineskip}
\end{figure}

Étant donné un choix de seuil~$\tau\in\R^+$ sur la persistance, on fusionne ensuite itérativement les clusters associés aux pics de persistance plus petite que~$\tau$ dans les clusters de leurs parents. Notons que, dans ce cas particulier, la structure stratifiée de l'espace~$G$ et la natu\-re linéaire par morceaux de la fonction~$\hat f$ garantissent que tout pic qui n'est pas un maximum global de~$\hat f$ sur la composante connexe de~$G$ où il se trouve a un parent. Les maxima globaux, quant à eux, ont par définition une persistance infinie. De ce fait, la procédure fusionne bien tous les clusters associés aux pics de persistance plus petite que~$\tau$ dans d'autres clusters, et \emph{in fine} dans des clusters asso\-ciés à des pics de persistance plus grande que~$\tau$. C'est ainsi que nous avons obtenu le clustering de la figure~\ref{fig:clustering} à partir de celui de la figure~\ref{fig:Koontz}, par un choix adéquat de seuil~$\tau$.

\subsubsection*{Choix du seuil de persistance}
En pratique, pour choisir la valeur du seuil~$\tau$ on peut s'appuyer sur le diagramme de persistance de~$\hat f$ dans~$G$, noté~$D(\hat f)$, que l'on a calculé lors du post-traitement. Grâce au théorème de stabilité~\ref{thm:stability_deg0}, on peut montrer que, sous des hypothèses d'échantillonnage adéquates, le diagramme~$D(\hat f)$ exhibe une séparation claire entre, d'une part, les intervalles de persistance des pics de~$\hat f$ correspondant aux pics de la densité sous-jacente~$f$, et d'autre part, les intervalles de persistance du reste des pics de~$\hat f$, qui correspondent à du bruit -- voir encore la figure~\ref{fig:ToMATo_PD} pour une illustration. Grâce à cette séparation, l'utilisateur (ou une méthode statistique) peut sélectionner un seuil~$\tau$ adapté. Les détails de l'analyse théorique et des garanties associées se trouvent dans l'article~\cite{cgos-pbc-13} qui a introduit cette approche, appelée \emph{ToMATo} pour \emph{Topological Mode Analysis Tool}.

\backmatter
\bibliographystyle{jepalpha+eid}
\bibliography{xups24-01}
\end{document}