TP Optimisation Différentiable

ENSAE, Avril 2018

Hicham Janati

import numpy as np
import scipy
from matplotlib import pyplot

Méthode de Newton

1.1 Introduction

Soit \(f\) une fonction de classe \(\mathcal{C}^1\) de \(\mathbb{R}^n\) dans \(\mathbb{R}^n\). Le but de la méthode de Newton est de résoudre \(f(x) = 0\). Soit \(x_0\) dans \(\mathbb{R}^n\). L’approximation de premier ordre de \(f\) dans un voisinage de \(x_0\) est donnée par: \[ \hat{f}(x) = f(x_0) + J(x_0)(x - x_0) \] Oû \(J\) est la matrice Jacobienne de f. Annuler la tangeante \(\hat{f}\) donne \(x = x_0 - J^{-1}f(x_0)\) et obtient donc la suite d’itérés: \[x_k = x_{k-1} - J^{-1}f(x_{k-1})\]

Question 1

Soit \(x^{\star}\) un zéro de \(f\). Supposons que \(J(x^{\star})\) est inversible et que \(f\) est de classe \(\mathcal{C}^2\). Montrez que la méthode de Newton a une convergence localement quadratique i.e qu’il existe une boule B centrée en \(x^{\star}\) telle que pour tout \(x_0\) dans B, il existe \(\alpha > 0\) telle que la suite de Newton vérifie: \[ \|x_{k+1} - x^{\star}\| < \alpha \|x_{k} - x^{\star}\|^2 \] ###### indication: Écrire l’approximation de deuxième ordre de f avec reste intégral.

---

Deux remarques importantes: - Newton est basée sur une approximation locale. La solution obtenue dépend donc du choix de \(x_0\). - \(J\) doit être inversible.

---

1.2 Optimisation sans contraintes

Soit \(g\) une fonction de classe \(\mathcal{C}^2\) de \(\mathbb{R}^n\) dans \(\mathbb{R}\).

Question 2

Adapter la méthode de Newton pour résoudre \(\min_{x \in \mathbb{R}^n} g(x)\).

Pour les questions 3-4-5, On prend \(g: x \mapsto \frac{1}{2}\|Ax - b\|^2 + \gamma \|x\|^2\), avec \(A \in \mathcal{M}_{m, n}(\mathbb{R})\), $b ^m $ et $ $

Question 3

Donner le gradient et la hessienne de g et complétez les fonctions gradient et hessian ci-dessous. Vérifiez votre gradient avec l’approximation numérique donnée par scipy.optimize.check_grad. ##### Question 4 Lorsque \(\gamma > 0\), montrez que la méthode de Newton converge en une itération indépendemment de \(x_0\).

Question 5

Complétez la fonction newton ci-dessous pour résoudre (2). Calculer l’inverse de la hessienne est très couteux (complexité \(O(n^3)\)), comment peut-on y remédier ? Vérifiez le point (4) numériquement.

seed = 1729 # Seed du générateur aléatoire
m, n = 50, 100
rnd = np.random.RandomState(seed) # générateur aléatoire
A = rnd.randn(m, n) # une matrice avec des entrées aléatoires gaussiennes
b = rnd.randn(m) # on génére b aléatoirement également 
gamma = 1.

def g(x):
    """Compute the objective function g at a given x in R^n."""
    Ax = A.dot(x)
    gx = 0.5 * np.linalg.norm(Ax - b) ** 2 + gamma * np.linalg.norm(x) ** 2
    return gx

def gradient_g(x):
    """Compute the gradient of g at a given x in R^n."""
    # A faire
    
def hessian_g(x):
    """Compute the hessian of g at a given x in R^n."""
    # A faire

Vous pouvez vérifier que votre gradient est bon en utilisant la fonction de scipy scipy.optimize.check_grad. Exécutez scipy.optimize.check_grad? pour obtenir la documentation de la fonction.

from scipy.optimize import check_grad
x_test = rnd.randn(n) # point où on veut évaluer le gradient
check_grad(g, gradient_g, x_test) # compare gradient_g à des accroissements petis de g

def newton(x0, g=g, gradient=gradient_g, hessian=hessian_g, maxiter=10, verbose=True):
    """Solve min g with newton method"""
    
    x = x0.copy()
    if verbose:
        strings = ["Iteration", "g(x_k)", "max|gradient(x_k)|"]
        strings = [s.center(13) for s in strings]
        strings = " | ".join(strings)
        print(strings)

    for i in range(maxiter):
        H = hessian(x)
        d = gradient(x)
        
        if verbose:
            obj = g(x)
            strings = [i, obj, abs(d).max()] # On affiche des trucs 
            strings = [str(s).center(13) for s in strings]
            strings = " | ".join(strings)
            print(strings)
        
        # A faire
        x = 

    return x

x0 = rnd.randn(n)
x = newton(x0)

1.3 Optimisation avec contraintes d’égalité

On s’intéresse à présent au problème avec contrainte linéaire: \[ \min_{\substack{x \in \mathbb{R}^n \\ Cx = d}} \frac{1}{2}\|Ax - b\|^2 + \gamma \|x \|^2 \]

Question 6

Donnez (en justifiant) le système KKT du problème.

Question 7

Expliquer comment peut-on utiliser la méthode de Newton pour résoudre le système KKT.

Question 8

Implémentez la fonction F dont on veut trouver un zéro et sa matrice Jacobienne.

Question 9

Implémentez la version de newton adaptée.

p = 5 # nombre de contraintes
C = rnd.randn(p, n)
d = rnd.randn(p)

def F(...):
    """Compute the function F."""
    # A faire
    
def jac_F(x):
    """Compute the jacobian of F."""
    # A faire

def newton_constrained( ):
    # A faire