import numpy as np
import scipy
from matplotlib import pyplot
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})$$
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 $$
Deux remarques importantes:
Soit $g$ une fonction de classe $\mathcal{C}^2$ de $\mathbb{R}^n$ dans $\mathbb{R}$.
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 \in \mathbb{R}^m $ et $\gamma \geq 0 $
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.
Lorsque $\gamma > 0$, montrez que la méthode de Newton converge en une itération indépendemment de $x_0$.
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)
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 $$
Donnez (en justifiant) le système KKT du problème.
Expliquer comment peut-on utiliser la méthode de Newton pour résoudre le système KKT.
Implémentez la fonction F dont on veut trouver un zéro et sa matrice Jacobienne.
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