TangenteX.com

Le problème

La recherche des extrema d'une fonction

La recherche des extrema, minimum ou maximum, d'une fonction d'une ou plusieurs variables est une activité très commune en physique numérique. S'agissant d'une fonction d'une variable, le problème est la plupart du temps simple, voire analytique. Mais lorsqu'on aborde les fonctions à plusieurs variables, ça se gâte ! Et comme généralement, les fonctions que l'on rencontre ont au moins deux variables...

Un rappel pour commencer. Soit f une fonction de \( \mathbb{R}^n \) dans \( \mathbb{R} \), continue et dérivable sur le domaine de travail. Maximiser la fonction f, c'est chercher le ou les points X de \( \mathbb{R}^n \) tels que \( \forall x \in \mathbb{R}^n \) alors \( f(X) > f(x) \). Bien sûr, chercher le minimum de la fonction f(x) revient à chercher le maximum de -f(x), ce qui fait que l'on peut se consacrer uniquement au problème de la recherche du maximum. Ou du minimum, mais j'ai choisi le maximum !

Dans cette page, je me contenterai de fonctions à deux variables. Mais les algorithmes que je vais mettre en oeuvre sont très facilement adaptables à une dimension quelconque.

Recherche du maximum global d'une fonction de deux variables

Pour illustrer cette page, j'ai donc choisi une fonction à deux variables dont la surface ne soit pas trop... ennuyeuse. Je vais utiliser une fonction que nous avons déjà rencontré ici : l'expression de l'intensité d'une onde plane diffractée par une ouverture carrée. Son expression est de la forme : \( \displaystyle \dfrac{I(x,y)}{I_0} = \left(\dfrac{sin(x)}{x}\right)^2 \left(\dfrac{sin(y)}{y}\right)^2\). Comme vous le constatez, j'ai laissé tomber les différents paramètres physiques qui ne m'intéressent pas ici.

Voici la surface représentative de cette fonction sur le domaine \( x \in [-\pi, \pi] \) et \( y \in [-\pi, \pi] \) :

Diffraction par une fente carrée

Cette surface comporte un maximum global sur ce domaine, que l'on peut calculer facilement et qui est situé au point X(0,0). Mais elle comporte aussi des maxima locaux, ce qui rend intéressant la recherche numérique du maximum global sur le domaine. On le voit mieux sur la représentation par courbes de niveau que voici :

Diffraction par une fente carrée

Ces deux figures sont tracées par le script Python Max2DTrace.py disponible ici.

Le but du jeu est donc de proposer une solution pour calculer et situer le maximum global de la fonction sans se trouver coincé sur un maximum local.

Il existe une solution analytique que vous connaissez sans doute : le calcul du gradient de la fonction f(x), s'il est possible. Pour rappel sur le gradient et sa définition, je vous renvoie à la page de TangenteX consacrée au sujet.

En effet, sur le domaine, le maximum global correspond au point X(x,y) tel que \(\displaystyle  \vec{grad} \: f(x,y) = 0 \), c'est à dire  \( \displaystyle  \vec{grad} \: f(x,y) = \dfrac{\partial f}{\partial x}(x,y)\vec{e}_x + \dfrac{\partial f}{\partial y}(x,y)\vec{e}_y = 0 \), soit encore \( \dfrac{\partial f}{\partial x}(x,y) = 0 \) ET \( \dfrac{\partial f}{\partial y}(x,y) = 0 \).

Il faut donc résoudre un système de deux équations différentielles aux dérivées partielles non linéaires. C'est faisable mais souvent très pénible. Aussi, les physiciens utilisent généralement des méthodes plus directes de recherche du maximum. Il en existe beaucoup de sortes, des algorithmes déterministes ou stochastiques. Je propose ici une méthode simple, basée sur l'utilisation du gradient, qui peut être facilement adaptée à des dimensions supérieures à 2. Ce n'est pas la plus précise, ni la plus efficace, mais son algorithme est facile à comprendre et les résultats suffisants pour nos besoins.

Une première tentative : la méthode du gradient simple

L'algorithme

Le principe est assez intuitif. On choisit un point appartenant au domaine. Puis on se déplace dans une direction à déterminer vers le point distant d'une grandeur qui reste aussi à déterminer. Et l'on se pose la question : ce point est-il plus "élevé" que le point d'où je viens ? Si oui, je continue dans la même direction et pour la même distance. Et je me repose la question. Si non, je change de direction et/ou de distance. Et je continue ainsi jusqu'à ce que je trouve le point le plus (le maximum de la fonction) ou bien que je sois fatigué !

L'utilisation de ce principe pour construire cet algorithme pose quelques questions : comment choisir le point initial ? Comment choisir la direction à suivre ? Comment déterminer la distance à parcourir entre deux points ? Est-on sûr d'avoir atteint le plus haut sommet ou bien s'est-on arrêté sur une colline ?

Je base mon algorithme sur une méthode simple pour déterminer la direction : la méthode du gradient, et je déterminerai, dans un premier temps, le point initial au lancement du script. Ce sera un paramètre que je pourrai faire varier pour évaluer son influence sur le résultat obtenu. Car pour répondre à la dernière question, il n'y aujourd'hui et à ma connaissance, aucun algorithme qui garantisse à 100% que l'on est pas bloqué sur un maximum local, du moins lorsqu'on parle de surfaces biscornues.

Comme vous le savez, le gradient d'une fonction f(x,y) en un point indique la pente de cette fonction dans les directions x et y. Bien sûr, je ne vais pas calculer analytiquement le gradient de f(x,y) mais je vais utiliser le développement de Taylor d'ordre 2 pour l'approximer numériquement. Cela me donne :

\( \displaystyle \dfrac{\partial f}{\partial x} = \dfrac{f(x + \delta,y) - f(x - \delta,y)}{2 \delta}\)

\( \displaystyle \dfrac{\partial f}{\partial y} = \dfrac{f(x, y + \delta) - f(x, y - \delta)}{2 \delta}\)

et donc, je peux aussi calculer sa norme :

\( \displaystyle \mid \vec{grad} \: f \mid = \sqrt{ \left( \dfrac{\partial f}{\partial x} \right)^2  + \left( \dfrac{\partial f}{\partial y} \right)^2 } \)

Moyennant ces deux informations, je calcule le déplacement unitaire que je vais appliquer pour me déplacer du point (xn,yn) ou point (xn+1,yn+1) :

\( \displaystyle x^{n+1} = x^n + \dfrac{\delta}{\mid \vec{grad} \: f \mid}  \dfrac{\partial f}{\partial x} \)

\( \displaystyle y^{n+1} = y^n + \dfrac{\delta}{\mid \vec{grad} \: f \mid}  \dfrac{\partial f}{\partial y} \)

Le paramètre \( \delta \) figure la distance parcourue, le pas, sur le plan (x,y) dans la direction du gradient.

Son implémentation en Python

L'implémentation en Python de l'algorithme que je viens de décrire ci-dessus est très simple. Je vous ferai grâce des appels aux packages de Python, ainsi qu'aux déclarations des variables et constantes, que vous retrouverez dans le script téléchageable ici. Je ne décrirai que les parties importantes du code.

La fonction à maximiser

Son équation est donnée par \( \displaystyle \dfrac{I(x,y)}{I_0} = \left(\dfrac{sin(x)}{x}\right)^2 \left(\dfrac{sin(y)}{y}\right)^2\). Vous remarquerez la forme \( \displaystyle \dfrac{sin(x)}{x} \), qui risque de vous poser problème pour le point (0,0). Votre calculateur n'aime pas les divisions par 0 ! Pourtant, le rapport \( \displaystyle \dfrac{sin(x)}{x} \) qui constitue ce que l'on appelle un sinus cardinal, est défini pour x = 0, en évaluant la limite du rapport \( \displaystyle \dfrac{sin(x)}{x} \) lorsque x tend vers 0. Par convention donc, sinc(0) = 1.

Aussi, deux possibilités s'offrent à nous pour implémenter notre équation :

J'utiliserai donc sinc().

def Fonction(x,y):
    z = fabs(sinc(x)*sinc(y))
    return z
Le code de l'algorithme
def MethodeGradientPF(x0,y0,pas,epsilon,maxiter):
    # initialisation des variables pour le calcul
    erreur = 2.*epsilon
    x = x0;y = y0;i = 0
    # traitement du point initial (0,0)
    if (x == 0.0) and (y == 0.0):
        X = 0.0; Y = 0.0; maximum = 1.0 
    # cas d'un point initial hors (0,0)
    else:       
        fmax = Fonction(x,y)
        # boucle de calcul du maximum. Les conditions de sortie sont :
        # soit la précision fixée par epsilon est atteinte,
        # soit le nombre maximum d'itérations est atteint.
        while ((erreur > epsilon) and (i < maxiter)):
            # calcul approché du gradient en (x,y)
            dx = Fonction(x + pas, y) - Fonction(x - pas, y)
            dy =
Fonction(x, y + pas) - Fonction(x, y - pas)
            gradF = hypot(dx,dy)
            # calcul des coordonnées du point (xn,yn) après déplacement
            x1 = x + pas*dx/gradF
            y1 = y + pas*dy/gradF
            # évaluation de la fonction au point (xn, yn)
            fmax1 = Fonction(x1,y1)
            # si f(xn,yn) est plus grand que f(x,y), alors le point est retenu comme un max possible
            if (fmax1 > fmax):
                x = x1;y = y1
                erreur = abs(fmax1 - fmax)
                fmax = fmax1
            # sinon, le pas est diminué
            else:
                pas /= 2.0
            i += 1
        # les conditions de sortie sont atteintes. On retourne les coordonnées du maximum et sa valeur.
        X = x; Y = y; maximum = fmax1
    return(X,Y,maximum,i)

Les commentaires dans le code détaillent les différentes étapes de l'algorithme.

Les résultats

Pour explorer les résultats obtenus avec l'algorithme décrit ci-dessus, je fixe la valeur des paramètres suivants :

epsilon = 1.e-6
maxiter = 50
pas = 0.1

Je me suis muni d'une impression de la courbe de niveau de la fonction sous test et j'ai lancé le script en fixant à chaque run un point initial différent dans le domaine \( x \in [-\pi, \pi] \) et \( y \in [-\pi, \pi] \) .

Commençons par vérifier le traitement du cas (0.0, 0.0). Comme attendu, le script retourne :

Maximum : 1.000  en (0.000,0.000) 0 iterations 

Choisissons maintenant un point dans le lobe principal de la courbe, par exemple en (0.5,0.5). Le script retourne :

Maximum : 1.000 en (-0.000,-0.000) 16 iterations

Le résultat est satisfaisant, mais il faut quand même 16 itérations pour l'obtenir...

Essayons maintenant en choisissant un point initial à proximité d'un lobe secondaire de la surface, par exemple en (0.7, 1.3), que j'ai repéré sur la courbe de niveau. J'obtiens, avec pas = 0.1 :

Maximum : 0.217 en (-0.000,1.431) 17 iterations

Mon algorithme s'est coincé sur le maximum local du lobe du haut ! Essayons, en partant du même point (0.7, 1.3) mais en utilisant un pas plus grand, dans mon cas pas = 1.0. J'obtiens maintenant :

Maximum : 1.000 en (-0.000,0.000) 19 iterations

L'algorithme trouve maintenant le bon résultat ! En fixant un pas plus grand, je permets un déplacement plus grand en début de recherche et donc j'évite de me retrouver coincé sur un maximum local.

Partons maintenant d'un point encore plus excentré, sur un lobe secondaire faible, par exemple en (0.0,2.4), toujours avec un pas de 1.0. J'obtiens :

Maximum : 1.000 en (0.000,0.000) 19 iterations

Là aussi, l'algorithme trouve le sommet principal et ne reste pas coincé sur le lobe secondaire. Si j'essaye avec un pas= 0.1, alors je reste sur le sommet du lobe secondaire, avec un maximum à 0.128 au point (0.000, 2.459).

Le résultat dépend donc de la géométrie de ma surface et du choix du pas de déplacement. Pas terrible, sauf si vous avez une idée de la "tête" de votre surface...

Une variante avec le choix aléatoire du pas et de la direction

Modification de l'algorithme

Nous venons de constater que le choix du point initial et la grandeur du pas de déplacement influaient fortement sur le résultat. Aussi, l'idée est de procéder à des essais multiples, en choisissant la direction et la grandeur du pas de déplacement de manière aléatoire, et en ne conservant que le plus grand maximum issu de ces essais. L'application de la méthode de Monte-Carlo à la recherche du maximum d'une fonction en quelque sorte...

Le pas de déplacement sera choisi selon une loi de distribution gaussienne centrée réduite. Nous avons vu plus haut l'influence de la taille du pas sur la possibilité de l'algorithme de rester coincé sur un maximum local. Avec une distribution gaussienne du pas, il arrivera qu'un grand pas soit généré, même avec une probabilité faible, ce qui décoincera notre algorithme. Le choix des caractéristiques de la gaussienne, moyenne et variance, doit répondre à deux contraintes : produire un pas de déplacement suffisamment petit pour que la précision de localisation du maximum soit de l'ordre de la précision attendue ; produire un pas suffisamment grand pour décoincer l'algorithme d'un maximum local. Ajoutons une autre contrainte : tous les déplacements qui provoqueraient une diminution de la valeur de la fonction ne seront pas pris en compte.

L'implémentation en Python

Seul la fonction implémentant l'algorithme est modfiée. Les commentaires dans le code devraient vous permettre de retrouver les instructions correspondantes aux modifications que je viens de décrire.

# méthode du gradient en partant d'un point initial (x0,y0) fixe,
# avec choix aléatoire du pas et de la direction du déplacement
def MethodeStochastique(x0,y0,pas,epsilon,maxiter):
    # initialisation
    erreur = 2.*epsilon
    x = x0;y = y0;
    i = 0
    # traitement du cas (0.0,0.)
    if (x == 0.0) and (y == 0.0):
        X = 0.0; Y = 0.0; maximum = 1.0 
    # cas général
    else:       
        fmax = Fonction(x,y)
        fmax1 = fmax/2.0
        # boucle principale de calcul       
        while ((erreur > epsilon) and (i < maxiter)):
            j = 0
            # boucle de choix aléatoire du pas et de la direction de déplacement
            # distribution gaussienne centrée réduite
            while (fmax1 <= fmax) and (j < nbtirages):     
                x1 = x + pas*gauss(mu, sigma)
                y1 = y + pas*gauss(mu, sigma)
                fmax1 = Fonction(x1,y1)
                j += 1
            if (fmax1 > fmax):
                x = x1;y = y1
                erreur = fabs(fmax1 - fmax)
                fmax = fmax1
            else:
                pas /= 2.0
            i += 1
        X = x; Y = y; maximum = fmax1   
    return(X,Y,maximum,i)

J'utilise la fonction gauss() du package random de Python pour générer aléatoirement les pas. Le script RechercheMax2DStoch.py est téléchargeable ici.

Les paramètres de l'algorithme sont :

# paramètres des tirages aléatoires   
nbtirages = 20
# loi normale centrée réduite
mu = 0.0        # moyenne de la distribution des pas
sigma = 1.0     # variance de la distribution des pas
# paramètres de l'algorithme
epsilon = 1.e-6
maxiter = 50
pas = 1.0

Les résultats

Commençons par vérifier le cas particulier (0.0,0.0) :

Maximum : 1.000  en (0.000,0.000) 0 iterations

Rien à signaler. Essayer maintenant un point proche d'un lobe secondaire (0.0,2.4) :

Maximum : 1.000  en (0.000,0.000) 24 iterations

et un autre (0.7, 1.3) :

Maximum : 1.000  en (-0.000,0.000) 23 iterations

Et maintenant un point à la limite du doamine (-3.1,3.1) :

Maximum : 1.000  en (0.000,0.000) 25 iterations

L'algorithme retourne des résultats corrects et précis en un nombre à peu près constant d'itérations. Parfois, nous sommes dans le domaine de l'aléatoire, le nombre d'itérations est supérieur (le max que j'ai constaté est de 30 itérations sur une vingtaine de runs) car l'algorithme se bloque sur un maximum local avant de repartir sur le tirage d'un pas un peu plus grand.

Contenu et design par Dominique Lefebvre - www.tangenteX.com octobre 2019 Licence Creative Commons Contact :