3.3. Méthode du recuit simulé

On suppose \(E\) fini et \(H:E \to \mathbf{R}\) une fonction réelle sur \(E\) à minimiser. Cette fonction \(H\) est appelée dans la suite fonction d’énergie. On note \(\underline H = \min_{x \in E} H(x)\) (l’énergie minimale) et \(\overline H = \max_{x \in E} H(x)\). L’algorithme de recuit simulé présenté ici est un algorithme stochastique (aléatoire) qui converge vers un minimum de \(H\) plus précisément un point de \(\operatorname{argmin} H = \big\{x \in E, H(x) = \underline H \big\}\) (il peut y avoir plusieurs éléments dans cet ensemble).

Une approche naïve consiste à parcourir tout l’espace \(E\) et de prendre un point de cet espace qui réalise un minimum \(\underline H\), mais si l’espace \(E\) est grand (très grand) cette approche n’est pas satisfaisante car trop coûteuse.

Une première idée est de considérer une loi de probabilité \(\mu_T\) (\(T\) un paramètre fixé) qui se concentre sur de points de \(E\) où \(H\) est proche de \(\underline H\) et d’utiliser l’algorithme de Metropolis pour construire une chaine de Markov sur \(E\) de loi invariante \(\mu_T\). La simulation de cette chaine de Markov est un algorithme qui parcourt \(E\) en se concentrant d’avantage sur des points de \(E\) où \(H\) est proche de \(\underline H\). Au fur et à mesure, on peut modifier la loi de probabilité pour qu’elle se concentre d’avantage autour de \(\operatorname{argmin}(H)\), ainsi on considère une famille \((\mu_{T_n})_{n \ge 1}\) de loi de probabilités où \((T_n)_{n \ge 1}\) est une suite qui tend vers 0. C’est le principe du recuit simulé.

3.3.1. Préliminaires: mesures de Gibbs

On présente d’abord la famille de loi de probabilité \((\mu_T)_{T > 0}\) qui se concentre sur \(\operatorname{argmin}(H)\) lorsque \(T\) tend vers 0.

Definition 3.6 (mesure de Gibbs)

La mesure de Gibbs associée à la fonction d’énergie \(H\) et à la température \(T > 0\) est la mesure de probabilité définie par

\[\begin{equation*} \forall x \in E, \quad \mu_T(x) = \frac{1}{Z_T} e^{-H(x) / T}, \end{equation*}\]

où \(Z_T = \sum_{x \in E} e^{-H(x) / T}\) est la constante de normalisation.

Cette section s’inspire très largement du chapitre 2 du livre [] dans lequel vous trouverez des compléments.

Dans la suite on note pour un sous-ensemble \(F\) de \(E\), \(\pi_F\) la loi de probabilité uniforme sur \(F\), c’est à dire

\[\begin{equation*} \forall x \in E, \quad \pi_F(x) = \frac{\mathbf{1}_F(x)}{\operatorname{card}(F)} \end{equation*}\]

Proposition 3.6

On a pour tout \(x \in E\),

\[\begin{equation*} \lim_{T \to 0} \mu_T(x) = \pi_{\operatorname{argmin}(H)}(x) \quad \text{et} \quad \lim_{T \to +\infty} \mu_T(x) = \pi_E(x). \end{equation*}\]

En particulier \(\lim_{T \to 0} \|\mu_T - \pi_{\operatorname{argmin}(H)}\| = 0\).

Proof. On prouve uniquement le cas \(T \to 0\) (le cas \(T \to \infty\) est similaire). Alors

\[\begin{equation*} e^{-H(x) / T} = e^{-(H(x) - \underline{H}) / T} e^{-\underline{H} / T} \mathbf{1}_{H(x) > \underline{H}} + e^{-\underline{H} / T} \mathbf{1}_{H(x) = \underline{H}}. \end{equation*}\]

Or comme \(E\) est fini, il existe \(\varepsilon > 0\) tel que \(\{H(x) > \underline H\} = \{H(x) \ge \underline H + \varepsilon\}\) donc sur cet ensemble on a \(e^{-(H(x) - \underline{H}) / T} \le e^{-\varepsilon / T}\) et

\[\begin{equation*} \lim_{T \to 0} e^{-(H(x) - \underline{H}) / T} e^{-\underline{H} / T} \mathbf{1}_{H(x) > \underline{H}} = 0. \end{equation*}\]

Ainsi

\[\begin{equation*} \lim_{T \to 0} e^{-H(x) / T} = \lim_{T \to 0} e^{-\underline H / T} \mathbf{1}_{\operatorname{argmin}(H)}(x), \end{equation*}\]

et en sommant sur tous les états

\[\begin{equation*} \lim_{T \to 0} Z_T = \sum_{x \in E} e^{-H(x) / T} = \lim_{T \to 0} e^{-\underline H / T} \operatorname{card}(\operatorname{argmin}(H)), \end{equation*}\]

d’où le résultat.

Proposition 3.7

On a pour tout \(T, T' > 0\)

\[\begin{equation*} \big\| \mu_T - \mu_{T'} \big\| \le \Big|\frac{1}{T} - \frac{1}{T'}\Big| (\overline H - \underline H). \end{equation*}\]

Proof. On suppose \(T \ge T'\) et \(\underline{H} = 0\) (quitte à poser \(\tilde H = H - \underline H\)). Ona

\[\begin{equation*} 0 \le e^{-H(x) / T} - e^{-H(x) / T'} = e^{-H(x) / T} \big(1 - e^{-(\frac{1}{T'} - \frac{1}{T}) H(x)} \big). \end{equation*}\]

En utilisant \(1 - e^{-x} \le x\) pour tout \(x \ge 0\) on prouve

\[\begin{equation*} e^{-H(x) / T} - e^{-H(x) / T'} \le \Big(\frac{1}{T'} - \frac{1}{T} \Big) \overline{H} e^{-H(x)/T}, \end{equation*}\]

d’où l’on déduit

\[\begin{equation*} |Z_T - Z_{T'}| \le \Big(\frac{1}{T'} - \frac{1}{T}\Big) \overline{H} Z_T, \end{equation*}\]

et en divisant par \({Z_T Z_{T'}}\)

\[\begin{equation*} \Big| \frac{1}{Z_T} - \frac{1}{Z_{T'}} \Big| \le \Big(\frac{1}{T'} - \frac{1}{T}\Big) \overline{H} \frac{1}{Z_{T'}}. \end{equation*}\]

Ainsi par la définition de la norme en variation totale

\[\begin{equation*} \big\|\mu_T - \mu_{T'}\big\| = \frac{1}{2} \sum_{x \in E} \Big|\frac{1}{Z_T} e^{-H(x) / T} - \frac{1}{Z_{T'}} e^{-H(x) / T'} \Big|. \end{equation*}\]

On introduit la quantité \(\frac{1}{Z_T} e^{-H(x) / T'}\) et par inégalité triangulaire on a

\[\begin{equation*} \big\|\mu_T - \mu_{T'}\big\| \le \frac{1}{2} \sum_{x \in E} \frac{1}{Z_T} \big| e^{-H(x) / T} - e^{-H(x)/T'}\big| + \frac{1}{2} \Big|\frac{1}{Z_T} - \frac{1}{Z_{T'}} \Big| \sum_{x \in E} e^{-H(x) / T'}. \end{equation*}\]

d’où l’on déduit \(\big\|\mu_T - \mu_{T'} \big\| \le \Big(\frac{1}{T'} - \frac{1}{T} \Big) \overline H\).

3.3.2. Algorithme

Soit \((T_n)_{n \ge 1}\) une suite déterministe positive qui décroît vers 0. Cette suite est appelée suite des températures et sera spécifiée dans la suite. Le choix de cette suite est extrêmement important en pratique et c’est un choix délicat. On considère \((\mu_n)_{n \ge 1}\) la suite des mesure de Gibbs associée à l’énergie \(H\) et à la température \(T_n\)

\[\begin{equation*} \mu_n = \mu_{T_n} = \frac{1}{Z_{T_n}} e^{-H(x) / T_n}. \end{equation*}\]

Par l’algorithme de Métropolis on construit une transition \(Q_n\) (ici il y a une dépendance en \(n\) mais cela ne change pas grand chose) qui est réversible par rapport à \(\mu_n\) et donc \(\mu_n\) invariante pour \(Q_n\). Soit \(P\) une matrice stochastique irréductible symétrique (qu’il faut choisir), alors

\[\begin{equation*} \mu_n Q_n = \mu_n \quad \text{avec} \quad Q_n(x, y) = \begin{cases} \Big( \frac{\mu_n(y)}{\mu_n(x)} \wedge 1 \Big) P(x, y) & \text{si $x \neq y$}, \\ 1 - \sum_{z \neq x} Q_n(x, z) & \text{sinon}. \end{cases} \end{equation*}\]

Par définition de \(\mu_n\) on a pour tous \(x \neq y\)

\[\begin{equation*} Q_n(x, y) = (e^{-(H(y) - H(x)) / T_n} \wedge 1) P(x, y), \end{equation*}\]

et on remarque que cette matrice dépend de \(H\) mais pas de \(Z_{T_n}\) la constante de normalisation qui est difficile à calculer car elle dépend de tout l’espace d’état \(E\). La matrice \(Q_n\) s’interprête de la façon suivante: on propose un voisin \(y\) de \(x\) avec la matrice \(P\)

• si \(H(y) < H(x)\) alors \((e^{-(H(y) - H(x)) / T_n} \wedge 1) = 1\) et on accepte toujours ce voisin: un voisin d’énergie plus basse est toujours accepté.

• si \(H(y) > H(x)\) alors on accepte ce voisin \(y\) d’énergie plus haute avec probabilité \(p_n = e^{-(H(y) - H(x)) / T_n}\): cette probabilité est d’autant plus petite que le saut d’énergie est élevé et que la température est basse mais cela permet à l’algorithme de mieux explorer l’espace \(E\) et d’éviter des minimas locaux. C’est une propriété essentielle de ce type d’algorithme.

La suite \((Q_n)_{n \ge 1}\) est une suite de matrice stochastique et on note \((X_n)_{n \ge 0}\) la chaine de Markov inohomogène de transition \((Q_n)_{n \ge 1}\) i.e.

\[\begin{equation*} \mathbf{P} \big(X_{n+1} = y | X_n = x\big) = Q_{n+1}(x, y). \end{equation*}\]

Si la loi initiale est \(X_0 \sim \nu_0\), on note \(\nu_n\) la loi de \(X_n\) et on vérifie aisément que \(\nu_n = \nu_0 Q_1 \dots Q_n\).

Intuitivement lorsque \(n\) tend vers l’infini, la loi de \(X_n\) est proche de \(\mu_n\) et \(\mu_n\) se concentre de plus en plus vers \(\pi_{\operatorname{argmin}(H)}\) la loi uniforme sur \(\operatorname{argmin}(H)\).

Proposition 3.8

Si \(P\) est symétrique irréductible et vérifie la condition de Doeblin et si

\[\begin{equation*} \forall n \ge 1, \quad T_n = \frac{\gamma}{\log(n)} \quad \text{avec} \quad \gamma > \max_{x,y} (H(y) - H(x)) \mathbf{1}_{P(x, y) > 0} \end{equation*}\]

alors pout toute loi initiaile \(\nu_0\) de \(X_0\) on a

\[\begin{equation*} \lim_{n \to +\infty} \big\| \nu_n - \mu_n \big\| = 0 \quad \text{où} \quad X_n \sim \nu_n. \end{equation*}\]

En particulier, \(\lim_{n \to +\infty} \mathbf{P}\big(H(X_n) > \underline H\big) = 0\).

Proof. On simplifie la preuve en supposant que \(P\) vérifie Doeblin avec \(l = 1\) c’est à dire qu’il existe \(\alpha > 0\) et \(c\) une probabilité tel que

\[\begin{equation*} \forall x, y \in E, \quad P(x, y) \ge \alpha c(y). \end{equation*}\]

On sait par Metropolis que \(\mu_n = \mu_n Q_n\) et \(\mu_{n+1} = \mu_{n+1} Q_{n+1}\) et par définition de \(\nu_n\) (la loi de \(X_n\)) que \(\nu_{n+1} = \nu_n Q_{n+1}\). Aini on a

\[\begin{equation*} \begin{aligned} \nu_{n+1} - \mu_{n+1} & = \nu_n Q_{n+1} - \mu_{n+1} Q_{n+1}, \\ & = (\nu_n - \mu_n) Q_{n+1} + (\mu_n - \mu_{n+1}) Q_{n+1}. \end{aligned} \end{equation*}\]

Ainsi en prenant la norme en variation totale on a

(3.10)\[ \big\|\nu_{n+1} - \mu_{n+1}\big\| \le \big\|(\nu_n - \mu_n) Q_{n+1}\big\| + \big\|(\mu_n - \mu_{n+1}) Q_{n+1}\big\|,\]

c’est à dire que l’erreur à l’itération \(n+1\) est majorée par la somme de 2 termes: le transport de l’erreur à l’itération \(n\) par le noyau \(Q_{n+1}\) de Metropolis et une erreur de biais entre 2 mesures de Gibbs.

Commençons par l’erreur de biais \(\big\|(\mu_n - \mu_{n+1}) Q_{n+1} \big\|\). Comme \(Q_{n+1}\) est une matrice stochastique on a

\[\begin{equation*} \big\|(\mu_n - \mu_{n+1}) Q_{n+1}\big\| \le \big\|\mu_n - \mu_{n+1}\big\| \end{equation*}\]

et d’après la proposition Proposition 3.7 on a

\[\begin{equation*} \big\|\mu_n - \mu_{n+1}\big\| \le \Big(\frac{1}{T_{n+1}} - \frac{1}{T_n}\Big) (\overline H - \underline H). \end{equation*}\]

Pour le terme \(\big\|(\nu_n - \mu_n) Q_{n+1}\big\|\). Soit \(\kappa = \max_{x, y} (H(y) - H(x)) \mathbf{1}_{P(x, y) > 0}\) qui vérifie notamment que \(a = e^{-\kappa / T_{n+1}}\) où \(a = \min_{x, y} \Big(\frac{\mu(y)}{\mu(x)} \mathbf{1}_{P(x, y) > 0} \Big)\). On vérifie que \(Q_{n+1}\) vérifie la condition de Doeblin avec \(l = 1\), la probabilité \(c\) et la constante \(\alpha_{n+1} = \alpha e^{-\kappa / T_{n+1}}\) et donc

\[\begin{equation*} \big\|(\nu_n - \mu_n) Q_{n+1}\big\| \le (1 - \alpha_{n+1}) \big\|\nu_n - \mu_n \big\|. \end{equation*}\]

Ainsi l’erreur à l’itération \(n+1\) vérifie d’après (3.10)

\[\begin{equation*} \big\|\nu_{n+1} - \mu_{n+1}\big\| \le (1 - \alpha_{n+1}) \big\|\nu_n - \mu_n\big\| + \beta_{n+1} \end{equation*}\]

en posant \(\beta_{n+1} = \Big(\frac{1}{T_{n+1}} - \frac{1}{T_n}\Big) (\overline H - \underline H)\). Pour obtenir un contrôle global à partir de ce contrôle local, on utilise le lemme de Gronwall discret que l’on rappelle dans cette preuve.

Lemma 3.1 (Gronwall discret)

Soit \((\alpha_n)_{n \ge 1}\) et \((\beta_n)_{n \ge 1}\) deux suites positives telles que \(\alpha_{n+1} \in ]0,1[\), \(\sum_n \alpha_n = +\infty\) et \(\lim_n \frac{\beta_n}{\alpha_n} = 0\). Si \((u_n)_{n \ge 0}\) est une suite positive vérifiant

\[\begin{equation*} \forall n \ge 0, \quad u_{n+1} \le (1 - \alpha_{n+1}) u_n + \beta_{n+1} \end{equation*}\]

alors \(\lim_{n \to +\infty} u_n = 0\).

On applique ce lemme avec \(\alpha_{n+1} = \alpha e^{-\kappa / T_{n+1}} = \alpha (n+1)^{- \kappa / \gamma}\) car \(T_n = \gamma / \log(n)\) et

\[\begin{equation*} \beta_{n+1} = \frac{1}{\gamma} (\log(n+1) - \log(n)) (\overline H - \underline H) = \frac{\overline H - \underline H}{\gamma} \log(1 + \frac{1}{n}). \end{equation*}\]

On vérifie que si \(\gamma > \kappa\) alors les hypothèses du lemme de Gronwall sont vérifiées et le résultat \(\lim_n \big\|\nu_n - \mu_n \big\| = 0\) est prouvé.

Proof. Gronwall discret On pose \(A_n = \frac{1}{\prod_{k=1}^n (1 - \alpha_k)}\) pour \(n \ge 1\) et \(A_0 = 1\). Alors

\[\begin{equation*} A_{n+1} u_{n+1} \le A_n u_n + A_{n+1} \beta_{n+1} \end{equation*}\]

Donc de l’écriture \(A_n u_n = A_0 u_0 + \sum_{k=1}^{n} (A_{k} u_k - A_{k-1} u_{k-1})\) on en déduit

(3.11)\[\begin{equation} u_n \le \frac{A_0 u_0}{A_n} + \frac{1}{A_n} \sum_{k=1}^n A_k \beta_k. \end{equation}\]

On vérifie que \(A_k = \frac{1}{\alpha_k} (A_k - A_{k-1})\) donc

(3.12)\[\begin{equation} u_n \le \frac{u_0}{A_n} + \frac{1}{A_n} \sum_{k=1}^n (A_k - A_{k-1}) \frac{\beta_k}{\alpha_k}. \end{equation}\]

Le premier terme de droite converge vers 0 car \(A_n\) tend vers l’infini et le second terme est une moyennisation à la Cesàro (pondérée par \(A_k - A_{k-1}\)) qui a même limite que la suite \((\frac{\beta_n}{\alpha_n})_{n \ge 1}\).