Algorithmes d'optimisation -- L3 MINT et doubles licences 2017/2018 -- Université Paris-Sud


TP 1: Minimisation de fonctionnelles quadratiques convexes

$\newcommand{\Rsp}{\mathbb{R}} \newcommand{\nr}[1]{\|#1\|} \newcommand{\abs}[1]{|#1|} \newcommand{\eps}{\varepsilon} \newcommand{\sca}[2]{\langle#1|#2\rangle} \newcommand{\D}{\mathrm{D}} \newcommand{\hdots}{\dots} \newcommand{\cond}{\mathrm{cond}}$ Dans ce TP, on s'intéresse à la minimisation sur $\Rsp^N$ de fonctionnelles de la forme suivante,

$$ \begin{aligned} f(x) &= \frac{1}{2} \sca{x}{Q x} + \sca{b}{x} \\ &= \frac{1}{2} \sum_{1\leq i\leq j \leq N} Q_{ij} x_i x_j + \sum_{1\leq i\leq N} b_i x_i, \end{aligned} $$

où $Q$ est une matrice symétrique définie positive et $b$ est un vecteur colonne. Notre objectif principal est de constater numériquement un phénomène expliqué en cours, à savoir que l'efficacité des méthodes de descente de gradient dépend crucialement du conditionnement de la matrice $Q$ (voir ci-dessous pour une définition).

Quelques commentaires sur les Notebook

Ce texte est rédigé sous la forme d'un notebook. Un notebook comporte des cellules de texte et des cellules de code, ici en Python. Quelques raccourcis clavier et remarques utiles:

  • CTRL+Entrée: exécute la cellule de code, et affiche son résultat.
  • Tab: Si l'on Tab après avoir tapé les premières lettres d'un nom de fonction, le système propose une liste de possibilités (ce qui peut permettre d'éviter des erreurs de frappe)
  • MAJ+Tab: Affiche la documentation sur la fonction. Très utile pour ne pas se tromper sur l'ordre des paramètres. On peut voir une documentation plus complète en cliquant sur le '+'.
  • CTRL+s: Enregistrer les modifications apportées au Notebook.
  • Le symbole [*] à côté d'une cellule de code indique que le noyau Python est toujours en train de calculer. On peut l'interrompre via Kernel -> Interrupt ou le redémarrer via Kernel -> Restart. Le noyau Python repart alors de zéro, et il faut donc relancer les cellules antérieures à celle sur laquelle on travaillait.

Une aide complète, ainsi que la documentation de Python et Numpy, est disponible dans le menu Aide.

Rappels de cours

Ce TP ne nécessite que quelques définitions et théorème du cours et des cours précédents, que l'on rappelle ici (le théorème de convergence ci-dessous sera démontré un peu plus tard, dans un cas plus général).

Proposition: Une fonction $f\in\mathcal{C}^2(\Rsp^d)$ est convexe si pour tout $x\in\Rsp^d$, $D^2 f(x)$ est une matrice symétrique positive, i.e. $\forall x\in\Rsp^d,\forall v\in \Rsp^d, \sca{v}{D^2 f(x) v} \geq 0$.

Proposition: Si $f \in \mathcal{C}^1(\Rsp^d)$ est convexe, alors $x^* = \arg\min_{x\in \Rsp^d} f(x) \Longleftrightarrow \nabla f(x^*) = 0. $

Définition: On appelle conditionnement d'une matrice symétrique définie positive $Q\in M_N(\Rsp)$ de valeurs propres $0< \lambda_1\leq \dots \leq \lambda_N$ la quantité $\cond(Q) = \lambda_N / \lambda_1$.

Théorème: Soit $f(x) = \frac{1}{2} x^T Q x + b^T x$ où $Q$ est une matrice symétrique définie positive, et soient $x^{(k)}$ la suite des itérées de l'algorithme de descente de gradient à pas optimal, c'est à dire: $$ \begin{cases} x^{(0)} \in \Rsp^d, \\ t^{(k)} = \arg\min_{t} f(x^{(k)} - t\nabla f_K(x^{(k)})) &\hbox{ pour $k\geq 0$}\\ x^{(k+1)} = x - t^{(k)} \nabla f_K(x^{(k)}). \end{cases} $$ alors, avec $x^* = \arg\min_x f(x)$, on a $$ f(x^{(k+1)}) - f(x^*) \leq c(f(x^{(k)}) - f(x^*)) \hbox{ où } c = 1 - \cond(Q)^{-1} < 1 $$

I. Gradient à pas optimal en dimension $N=2$

Pour commencer, on commence par considérer la minimisation de la fonction $f_K:\Rsp^2\to\Rsp$ définie par $f_K(x) = \frac{1}{2}Kx_1^2 + \frac{1}{2}x_2^2$ où $K$ est une constante strictement positive.

QI.1. Calculer le gradient et la matrice hessienne de $f_K$, puis montrer que $f_K$ est convexe. On note $x^* = (0,0)$ le minimiseur global de $f_K$.

QI.2. Étant donné un point $x = (x_1,x_2)$ et $d = -\nabla f_K(x)$, calculer le pas optimal, i.e. l'unique minimiseur de la fonction $g: t \mapsto f_K(x + t d)$.

QI.3. Programmer l'algorithme du gradient à pas optimal, et lel tester pour $K=2.$: $$ \begin{cases} x^{(0)} = (1,K) \\ x^{(k+1)} = x - t(x^{(k)}) \nabla f_K(x^{(k)}). \end{cases} $$ On arrêtera les itérations lorsque $f_K(x^{(k)}) \leq f_K(x^*) + 10^{-8}$. On stockera la liste des itérées $x^{(k)}$ et des valeurs $f(x^(k))$ comme suggéré ci-dessous. Enfin, on tracera ensuite sur deux figures : la trajectoire des itérées $(x_1^{(k)}, x_2^{(k)})$ (via la fonction plt.plot) et l'évolution $(k, \log(f(x^{(k)}) - f(x^*)))$ (on utilisera la fonction plt.semilogy). Que se passe-t-il lorsqu'on change $K$ ?

In [5]:
# on importe les modules numpy et pyplot
import numpy as np
import matplotlib.pyplot as plt
# les deux commandes suivante paramètrent l'affichage des figures
%matplotlib inline
plt.rcParams['figure.figsize'] = [9.,6.]

# initialiser les variables K,Q, x0, fmin
K = 2.
X = [] # tableau pour stocker la liste des itérées x^k
F = [] # tableau pour stocker la liste des valeurs f(x^k)
x = np.array([1.,K]) 

while(True):
    # calculer la direction de descente d, le pas t et arrêter si le critère d'arrêt est vérifié
    # <completer>
    F.append(f)
    X.append(x)
    x = x + t*d
X = np.array(X)
F = np.array(F)

# afficher les résultats
# <completer>

QI.4. Refaire l'expérience en choisissant $K = 10,100,500$. Calculer le taux de décroissance moyen (c'est-à-dire la moyenne de $(f_K(x^{(k+1)})-f(x^*))/(f_K(x^{(k)}) - f(x^*)$) et le comparer à la borne donnée dans le théorème rappelé au début du sujet.

In [5]:
# <completer>
        # <completer>
        x = x + t*d
        F.append(f)
    F = np.array(F)
    niter = len(F)
    print('nombre d\'itérations pour K=%d: %d' % (K, len(F)))
    print('taux de décroissance moyen pour K=%d: %g' % (K,np.mean(F[1:niter]/F[0:(niter-1)])))
    print('taux de décroissance théorique K=%d: %g\n' % (K,1-min(1.,K)/max(1.,K)))

II. Gradient à pas optimal pour un problème de moindres carrés

On commence par écrire l'algorithme de gradient à pas optimal pour minimiser une fonction quadratique convexe sur
$\Rsp^N,$ de la forme $$f(x) = \frac{1}{2} \sca{Q x}{x} + b^T x $$ où $Q$ est une matrice $N\times N$ symétrique définie positive et $b\in\Rsp^N$ est un vecteur.

QII.1 Calculer le gradient puis la Hessienne de $f$. En déduire que $f$ est convexe, puis l'équation vérifiée par le minimiseur $x^*$. Cette équation nous sera utile pour vérifier le fonctionnement de nos programmes.

QII.2 Étant donné $x\in \Rsp^N$ et $d = - \nabla f(x)$, calculer le pas optimal $t(x)$, i.e. le minimiseur de $g: t\mapsto f(x + td)$.

QII.3 Programmer une fonction gradient_optimal(Q,b,x0,err), prenant en argument la matrice $Q$, le vecteur $b$ et le point de départ $x^{(0)}$, qui calculera les itérées

$$\begin{cases} d^{(k)} = - \nabla f(x^{(k)}) \\ x^{(k+1)} = x^{(k)} + t(x^{(k)}) d^{(k)} \end{cases} $$

Les itérations seront interrompues lorsque $\nr{d^{(k)}} \leq $ err (on pourra utiliser la fonction np.linalg.norm pour calculer la norme). Elle retournera le dernier point $x^{(k)}$ trouvé ainsi que deux vecteurs $E,F$ tels que $E^{(k)} = \nr{d^{(k)}}$ et $F^{(k)} = f(x^{(k)})$. Tester cette fonction avec une matrice $Q = A^T A + \mathrm{Id}$, où $A$ est une matrice aléatoire et $b$ est un vecteur aléatoire, pour $N=10$. On comparera la solution construite par l'algorithme de descente de gradient à la solution exacte $-Q^{-1}b retournée par -np.linalg.solve(Q,b).

In [6]:
def gradient_optimal(Q,b,x0,err=1e-6):
    x = x0
    niter=0
    E = []
    F = []
    
    k = 0 # nombre d'itérations
    while (True): 
        k = k+1
        if k > 1e6: # maximum de 10^6 itérations
            print('erreur: nombre maximum d\'itérations atteint')
            break
        # calculer la direction de descente
        # <completer>
        # vérifier le critère d'arrêt, et quitter la boucle (avec break) s'il est vérifié
        # <completer>
        # calculer le pas de descente et mettre à jour x
        # <completer>
    E = np.array(E)
    F = np.array(F)
    return x,E,F
In [7]:
# test de la fonction gradient_optimal pour un Q,b aléatoire
N = 10
A = np.random.randn(N,N)
Q = np.dot(A.T,A)+np.eye(N,N)
b = np.random.rand(N)
# on calcule la solution exacte en résolvant le système linéaire Qx = -b (cf question précédente)
sol = -np.linalg.solve(Q,b)
# appliquer le gradient à pas optimal et vérifier qu'on obtient la même solution
# <completer>

On s'intéresse maintenant à un problème de moindres carrés: on a un jeu de donnée constitué de points $(x_i,y_i)_{1\leq i\leq n} \in\Rsp^2$. On cherche à trouver $(\mu_0, \mu_1) \in\Rsp^2$ minimisant la fonction

$$ f(\mu) = \frac{1}{2} \sum_{1 \leq i \leq n} |\mu_0 + \mu_1 x_i - y_i|^2. $$

Ce problème d'optimisation peut être compris de la manière suivante: étant donnée une abscisse à l'origine $\mu_0$ et une pente $\mu_1$, on peut calculer $\eps_i = |\mu_0 + \mu_1 x_i - y_i|$, qui correspond à la distance euclidienne entre le la donnée $(x_i,y_i)$ et le point correspondant sur la droite $(x_i,\mu_0 + \mu_1 x_i)$. Le problème d'optimisation cherche donc à trouver l'abscisse à l'origine $\mu_0$ et la pente $\mu_1$ qui minimise la somme des erreurs $\eps_1^2 + \hdots + \eps_n^2$, et assure donc que la droite correspondante passe au plus près du jeu de données. Cette méthode a été inventée simultanément par Gauss et Legendre.

Pour nos expériences, les points de données $(x_i,y_i)$ sont sur le graphe de la fonction cosinus, perturbée par l'ajout d'un bruit. Dans la suite, on choisira $x_i = -1+2\frac{i-1}{n-1}$ où $1\leq i\leq n$.

In [8]:
n=50
X = np.linspace(-1,1,n);
Y = np.sin(np.pi*X) + .1*np.random.rand(n)
plt.plot(X,Y,'.',label='donnees')
plt.legend()
Out[8]:
<matplotlib.legend.Legend at 0x10c32dc18>

QII.4 Construire une matrice $A$ ayant $n$ lignes et $2$ colonnes telle que le vecteur $Z := (\mu_0 + \mu_1 x_i)_{1\leq i\leq n}$ s'écrive $Z = A\mu$ (produit matrice vecteur). Remarquer que l'on peut alors écrire $f(\mu) = \frac{1}{2} \nr{A \mu - Y}^2$ où $Y = (y_i)_{1\leq i\leq n}$.

QII.5 Déterminer $Q \in M_2(\Rsp),$ $b\in\Rsp^2$ et $c\in\Rsp$ tels que $f(\mu) = \frac{1}{2} \sca{\mu}{Q \mu} + \sca{b}{\mu} + c$. Démontrer que $A$ est injective, puis en déduire que la matrice $Q$ est symétrique définie positive. En déduire la solution exacte $\mu^*$ du problème de minimisation de $f$ vérifie $A^T A \mu = A^T \mu.

QII.6 Calculer le minimum $\mu^*$ de $f$ à l'aide de la fonction gradient_optimal, tracer sur une même figure les points $(X,Y)$ et la droite $(X, A\mu^*)$. Vérifier la correction du résultat en le comparant à celui obtenu par résolution d'un système linéaire (cf question II.1).

In [12]:
# calculer A, Q, b, puis appliquer l'algorithme du gradient à pas optimal
# <completer>

QII.7 On s'intéresse maintenant à une généralisation du problème de moindres carrés, où l'on cherche à approcher le jeu de données $(x_i,y_i)$ par $(x_i,P(x_i))$ où $P$ est un polynôme. On considère le problème d'optimisation: $$ \min_{\mu_0,\hdots,\mu_d} f_d(\mu) := \frac{1}{2} \nr{A_d \mu - Y}^2 \quad \hbox{ où } A_d = \begin{pmatrix} 1 & x_1 & x_1^2 & \hdots & x_1^d \\ 1 & x_2 & x_2^2 & \hdots & x_2^d \\ \vdots & \vdots & \vdots & &\vdots \\ 1 & x_n & x_n^2 & \hdots & x_n^d \end{pmatrix}$$ Montrer que ce problème est équivalent au problème suivant, et l'interpréter

$$ \min_{P \in \Rsp_d[X]} \frac{1}{2} \sum_{1\leq i\leq N} \nr{P(x_i) - y_i}^2, $$

QII.8 Mettre le problème de minimisation sous la forme $f_d(\mu) = \frac{1}{2}\mu^T Q_d \mu + b_d^T \mu + c_d$ où $Q_d,b_d,c_d$ sont à déterminer. Démontrer que la matrice $A_d$ est injective lorsque $d\leq n-1$ (où $n$ est le nombre de points), puis en déduire que $Q$ est symétrique définie positive.

QII.9 Résoudre le problème d'optimisation pour $2\leq d \leq 7$ via la fonction gradient_optimal, en fixant err=1e-4. Interpréter l'accroissement du nombre d'itérations de l'algorithme en calculant le conditionnement des matrices $A_d$ via la fonction np.linalg.cond.

In [13]:
# on commence par construire la matrice A pour d=2, en guise d'exemple:
A = np.vstack([X**0, X**1, X**2]).T
# ou, de manière équivalente (il faut passer du temps à comprendre cette construction...)
A = np.vstack([X**i for i in range(3)]).T

# construire la matrice Q2, puis résoudre le problème d'optimisation par gradient à pas optimal
# <completer>
# on applique ensuite l'algorithme du gradient à pas optimal pour retrouver mu
mu,E,F = gradient_optimal(Q,b,np.zeros(3),1e-4)


# une fois que ça fonctionne pour d=2, on traite les autres cas dans une boucle
for d in range(2,8):
    # calculer les matrice A et Q, le vecteur b
    # <completer>
    # on applique ensuite l'algorithme du gradient à pas optimal pour retrouver mu
    mu,E,F = gradient_optimal(Q,b,np.zeros(d+1),1e-4)
    print('deg=%d, condition=%f,niter=%d' % (d,np.linalg.cond(Q),len(E)))
    plt.subplot(2,3,d-1) # on organise les figures dans un tableau à 2 lignes et 3 colonnes
    plt.plot(X,np.dot(A,mu))
    plt.plot(X,Y,'.')

QII.10 On pose $\tilde{Q}_d = Q_d + \gamma \mathrm{Id}$, où $Q_d$ est la matrice de la question précédente et $\gamma$ est une constante positive (par exemple $\gamma=0.1$ ou $0.01$), et on pose $\tilde{f}_d(\mu) := \frac{1}{2} \mu^T \tilde{Q}_d \mu + b_d^T \mu + c_d$. Démontrer que $\cond(\tilde{Q}_d) < \cond(Q_d)$, et vérifier que l'algorithme du gradient à pas optimal converge en moins d'itérations.

In [20]:
# NB: ci-dessous, on peut utiliser la fonction np.eye() pour construire la matrice identité
# <completer>

III Exercice (valeurs propres et conditionnement)

QIII.1 Soit $Q$ une matrice symétrique de taille $d\times d$ et $\lambda_1\leq \hdots\leq \lambda_d$ ses valeurs propres. Démontrer que pour tout $x \in \Rsp^d$, $$ \lambda_1 \nr{x}^2 \leq \sca{x}{Qx} \leq \lambda_d \nr{x}^2. $$ En déduire que $$ \lambda_1 = \min_{x \neq 0} \frac{\sca{x}{Qx}}{\nr{x}^2} \qquad \lambda_d = \max_{x \neq 0} \frac{\sca{x}{Qx}}{\nr{x}^2} $$

QIII.2 On considère maintenant la matrice $Q_d = A_d^T A_d$ où $A_d$ est définie en QII.7. En considérant $\mu = (1,0,\hdots,0)$ (resp. $\mu = (0,\hdots,0,1)$) démontrer que $$\lambda_d \geq n \qquad (\hbox{resp. } \lambda_1 \leq \sum_{1\leq i\leq n} x_i^{2d})$$ En déduire une minoration de $\cond(Q_d)$ ne faisant intervenir que $n$ et $d$.