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$.

Réponse: Le calcul des dérivées partielles donne $$ \nabla f_K(x) = \begin{pmatrix} K x_1\\ x_2\end{pmatrix} \qquad \D^2 f_K(x) = \begin{pmatrix} K & 0 \\ 0 & 1 \end{pmatrix}$$ Comme la Hessienne est positive, la fonction est convexe (elle est même strictement convexe comme $\D^2 f_K(x)$ est définie positive).

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)$.

Réponse: Il suffit de calculer $g$ et sa dérivée et de trouver l'unique $t$ tel que $g'(t) = 0$: $$ g(t) = \frac{1}{2} K (x_1 - t K x_1)^2 + \frac{1}{2} (x_2 - tx_2)^2 $$ $$ g'(t) = K^2 x_1 (t K x_1 - x_1) + x_2(tx_2 - x_2) = t(K^3 x_1^2 + x_2^2) - (K^2 x_1^2 + x_2^2) $$ soit $$ t^* = (K^2 x_1^2 + x_2^2)/(K^3 x_1^2 + x_2^2) $$

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é
    d = -np.array([K*x[0], x[1]]) 
    f = .5*(K*x[0]**2 + x[1]**2)  
    if f <= 1e-8:
        break
    t = ((K**2)*x[0]**2 + x[1]**2)/((K**3)*x[0]**2 + x[1]**2)
    F.append(f)
    X.append(x)
    x = x + t*d
X = np.array(X)
F = np.array(F)

# afficher les résultats
plt.plot(X[:,0],X[:,1])
plt.axis('equal')
plt.figure()
plt.semilogy(range(len(F)),F)
Out[5]:
[<matplotlib.lines.Line2D at 0x10b11f128>]

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]:
for K in [10,100,500]:
    F = [] # tableau pour stocker la liste des valeurs f(x^k)
    x = np.array([1.,K]) 

    for i in range(5000): # on limite à 5000 itérations
        # calculer la direction de descente d, le pas t et arrêter si le critère d'arrêt est vérifié
        d = -np.array([K*x[0], x[1]]) 
        f = .5*(K*x[0]**2 + x[1]**2)  
        if f <= 1e-8:
            break
        t = ((K**2)*x[0]**2 + x[1]**2)/((K**3)*x[0]**2 + x[1]**2)
        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)))
nombre d'itérations pour K=10: 56
taux de décroissance moyen pour K=10: 0.669421
taux de décroissance théorique K=10: 0.9

nombre d'itérations pour K=100: 674
taux de décroissance moyen pour K=100: 0.960788
taux de décroissance théorique K=100: 0.99

nombre d'itérations pour K=500: 3770
taux de décroissance moyen pour K=500: 0.992032
taux de décroissance théorique K=500: 0.998

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.

Réponse: On développe $f(x+v)$ pour obtenir $$ f(x+v) - f(x)= \sca{Q x}{v} + \sca{b}{v} + \frac{1}{2}\sca{Qh}{h} + o(\nr{v}^2), $$ ce qui donne par définition du gradient $\nabla f(x) = Qx + b$. Par identification, on trouve également $D^2 f(x) = Q$. Ces deux résultats peuvent également être trouvés par un calcul de dérivées partielles. Tout minimiseur de $f$, $x^*$ vérifie l'équation $\nabla f(x^*) = 0$. Cette équation n'a qu'une seule solution, $x^* = -Q^{-1}b$.

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)$.

Réponse: Par dérivation de fonctions composées, on a $$g'(t) = \sca{d}{\nabla f(x+td)} = \sca{d}{Q(x+td) + b}.$$ Ainsi, $g'(t) = 0$ si $t\sca{d}{Q d} = -\sca{d}{Qx + b} = \sca{d}{d}$, i.e. $$ t = \frac{\sca{d}{d}}{\sca{d}{Qd}}. $$

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
        d = -(np.dot(Q,x)+b)
        E.append(np.linalg.norm(d))
        F.append(.5*np.dot(x,np.dot(Q,x)) + np.dot(b,x))
        # vérifier le critère d'arrêt, et quitter la boucle (avec break) s'il est vérifié
        if np.linalg.norm(d) <= err:
            break
        # calculer le pas de descente et mettre à jour x
        t = np.linalg.norm(d)**2/np.dot(np.dot(Q,d), d)
        x = x + t*d
    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
x,E,F = gradient_optimal(Q,b,np.zeros(N),1e-8)
print(np.linalg.norm(x-sol))
7.87036406235e-09

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}$.

Réponse: On pose $A = [1 X]$ la matrice dont la première colonne est constituée entièrement de $1$ et la deuxième est le vecteur $X = (x_1,\hdots,x_n)$. Alors, $Z_i = \mu_0 + \mu_1 x_i = (A\mu)_i$.

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.

Réponse: En développant, on a $$ \begin{aligned} f(\mu) &= \frac{1}{2} \nr{A \mu - Y}^2 \\ &= \frac{1}{2} (A\mu-Y)^T (A\mu - Y) \\ &= \frac{1}{2}\mu^T (A^T A) \mu - \frac{1}{2}(A\mu)^T Y - \frac{1}{2}Y^T (A\mu) + \frac{1}{2}Y^T Y\\ &= \frac{1}{2}\mu^T Q \mu + b^T \mu + c \end{aligned} $$ où $Q = A^T A$, $b = - A^T Y$ et $c = \frac{1}{2}Y^T Y$.

Pour montrer que $A$ est injective, il suffit de vérifier que n'importe quelle matrice 2x2 extraite est inversible. Montrons maintenant que $Q$ est symétrique définie positive:

$$ \mu^T Q \mu = \mu^T A^T A \mu = (A \mu)^T A\mu = \nr{A\mu}^2. $$

Ainsi, $\mu^T Q \mu \geq 0$ avec égalité si et seulement si $A\mu = 0$, soit $x = 0$. Finalement, on sait que la solution $\mu^*$ vérifie $\nabla f(\mu^*) = 0$, soit $Q\mu^* = -b = A^T Y$.

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
A = np.vstack([np.ones(n),X]).T
Q = np.dot(A.T,A)
b = -np.dot(A.T,Y)
mu,E,F = gradient_optimal(Q,b,np.zeros(2))
plt.plot(X,Y,'.')
plt.plot(X,np.dot(A,mu))
Out[12]:
array([[  5.00000000e+01,  -3.77475828e-15],
       [ -3.77475828e-15,   1.73469388e+01]])

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, $$

Réponse: On identifie $\mu$ comme le vecteur des coefficients du polynôme $P$ : $P(X)=\sum^d_{i=0} a_iX^i$. On en déduit que $A_d\mu$ est le vecteur colonne ($n$ lignes) des évaluations de $P$ en $x_i$ et donc $\sum_{1\leq i\leq n}\nr{P(x_i) - y_i}^2= \nr{A_d \mu - Y}^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.

Réponse: Il n'y a pas de différence avec le cas $d=1$ : En développant, on a $$ \begin{aligned} \frac{1}{2}f(\mu) &= \frac{1}{2}\nr{A_d \mu - Y} \\ &= \frac{1}{2}(A_d\mu-Y)^T (A_d\mu - Y) \\ &= \frac{1}{2}\mu^T (A_d^T A_d) \mu - \frac{1}{2}(A_d\mu)^T Y - \frac{1}{2}Y^T (A_d\mu) + \frac{1}{2}Y^T Y\\ &= \frac{1}{2}\mu^T Q_d \mu + b_d^T \mu + c_d \end{aligned} $$ où $Q = A^T A$, $b = - A^T Y$ et $c = \frac{1}{2}Y^T Y$.

Vu les dimensions, l'opérateur $A_d$ ne peut pas être injectif si $d\geq n$. Soit $\mu$ tel que $A_d \mu = 0$, et $P(X) = \sum_{0\leq i\leq d} \mu_i X^i$. Alors $P(x_i) = 0$ pour tout $1\leq i\leq n$. Comme tout les $x_i$ sont distincts et que $n\geq d+1$, cela implique que $P = 0$, soit $\mu = 0$. Montrons maintenant que $Q$ est symétrique définie positive:

$$ x^T Q_d x = x^T A_d^T A_d x = (A_d x)^T A_d x = \nr{A_d x}^2. $$

Ainsi, $x^T Q_d x \geq 0$ avec égalité si et seulement si $A_d x = 0$, soit $x = 0$ par injectivité.

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
Q = np.dot(A.T,A)
b = -np.dot(A.T,Y)
# 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
    A = np.vstack([X**i for i in range(0,d+1)]).T
    Q = np.dot(A.T,A)
    b = -np.dot(A.T,Y)
    # 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,'.')
deg=2, condition=13.312869,niter=41
deg=3, condition=62.230914,niter=329
deg=4, condition=325.519883,niter=371
deg=5, condition=1666.465961,niter=5936
deg=6, condition=9034.326437,niter=10802
deg=7, condition=48152.315204,niter=73234

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.

Réponse: Soient $\lambda_1 \leq \hdots \leq \lambda_d$ les valeurs propres de $Q_d$. Alors, les valeurs propres de $\tilde{Q}_d$ sont $\lambda_1+\gamma\leq\hdots\leq\lambda_d+\gamma$. On a donc $$ \cond(Q_d) =\frac{\lambda_d}{\lambda_1} \qquad \cond(Q_d) =\frac{\lambda_d+\gamma}{\lambda_1+\gamma}. $$ Il suffit donc de démontrer que si $\gamma\geq 0$ et $\lambda_i\geq 0$, $$ \frac{\lambda_d+\gamma}{\lambda_1+\gamma} \leq \frac{\lambda_d}{\lambda_1}, $$ ce qu'on vérifie facilement.

In [20]:
# NB: ci-dessous, on peut utiliser la fonction np.eye() pour construire la matrice identité
plt.plot(X,Y,'.')
for d in range(2,8):
    A = np.vstack([X**i for i in range(0,d+1)]).T
    Q = np.dot(A.T,A) + 0.01*np.eye(d+1)
    b = -np.dot(A.T,Y)
    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,'.')
deg=2, condition=13.283965,niter=33
deg=3, condition=61.564747,niter=334
deg=4, condition=308.754186,niter=619
deg=5, condition=1303.272224,niter=4210
deg=6, condition=3672.104974,niter=4518
deg=7, condition=5481.963218,niter=9726

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} $$

Réponse: Comme $Q$ est symétrique, elle peut être diagonalisée en base orthonormale. Il existe donc $P$ une matrice orthogonale (telle que $P^T P = I_d$) et $D = \mathrm{diag}(\lambda_1,\hdots,\lambda_d)$ une matrice diagonale telle que $Q = P^T D P$. Alors, $$\sca{x}{Qx} = \sca{x}{P^T D P x} = \sca{P x}{D P x}. $$ On pose $y = P x = (y_1,\hdots, y_d)$. Alors, comme $\lambda_d\geq \lambda_i$, $$ \sca{P x}{D P x} = \sum_{1 \leq i\leq d} \lambda_i y_i^2 \leq \lambda_d \sum_{1\leq i\leq d} y_i^2 = \lambda_d \nr{Px}^2 $$ Or, comme $P$ est orthogonale, $\nr{Px}^2 = \sca{Px}{Px} = \sca{P^T P x}{x} = \nr{x}^2$. On en déduit l'inégalité $$ \sca{x}{Qx}\leq \lambda_d \nr{x}^2, $$ soit $$ \max_{x \neq 0} \frac{\sca{x}{Qx}}{\nr{x}^2} \leq \lambda_d. $$ On vérifie qu'il y a bien égalité en prenant $x$ un vecteur propre de $Q$ associé à la valeur propre $\lambda_d$. Le cas $\lambda_1$ se traite exactement de la même manière.

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$.

Réponse: On applique la deuxième inégalité de la question précédente, en prenant $\mu=(1,0,\hdots,0)$:

$$ \lambda_d \geq \frac{\sca{Q_d \mu}{\mu}}{\nr{\mu}^2} = \sca{A_d^T A_d \mu}{\mu} = \nr{A_d\mu}^2 = n, $$

et le même type de raisonnement montre que $\lambda_1 \leq \sum_{1\leq i\leq n} x_i^{2d}$. De plus,

$$ \begin{aligned} \sum_{1\leq i\leq n} x_i^{2d} &= 2n \left(\frac{1}{2n} \sum_{1\leq i\leq n} x_i^{2d}\right) &\sim_{n\to+\infty} 2n \int_{-1}^1 x^{2d} d x = \frac{4n}{2d + 1}, \end{aligned} $$

où l'on a utilisé la convergence des sommes de Riemann. Il faut noter qu'on pourrait être plus précis en découpant l'intégrale entre $[-1,0]$ et $[0,1]$, où la fonction $x\mapsto x^{2d}$ est monotone. On obtient donc

$$ \cond(Q_d) \geq \frac{\lambda_d}{\lambda_1} \sim_{n\to +\infty} \frac{2d+1}{4}. $$

Ainsi, le conditionnement se déteriore lorsque le degré croît. On remarquera que la borne donnée ci-dessus n'est pas très bonne, en comparant avec les calculs numériques.