import plotly
import numpy as np
import plotly.plotly as py
import plotly.graph_objs as go
import matplotlib.pyplot as plt
import scipy.io as sio
from scipy import stats
%matplotlib inline
plotly.offline.init_notebook_mode(connected=True)
from IPython.core.display import display, HTML
# The polling here is to ensure that plotly.js has already been loaded before
# setting display alignment in order to avoid a race condition.
display(HTML(
'<script>'
'var waitForPlotly = setInterval( function() {'
'if( typeof(window.Plotly) !== "undefined" ){'
'MathJax.Hub.Config({ SVG: { font: "STIX-Web" }, displayAlign: "center" });'
'MathJax.Hub.Queue(["setRenderer", MathJax.Hub, "SVG"]);'
'clearInterval(waitForPlotly);'
'}}, 250 );'
'</script>'
))
Dans ce TP on va chercher à implémenter certaines notions vues en cours concernant l'optimisation convexe. Le premier exercice a pour but de montrer que la vitesse de convergence d'un algorithme est considérablement influencée par certaines caractéristiques du problème. On implémentera ensuite une méthode de Newton dans le deuxième exercice. Le dernier exercice cherche à montrer que les conditions de KKT permettent souvent de donner une interprétation aux variables duales.
Remarque préliminaire : Les vitesses de convergence des algorithmes de se représentent généralement sur un graphe logarithmique, il faut donc penser à utiliser les fonctions semilogy, loglog de Python
Le but de ce problème est d'implémenter une méthode de Newton dans le cas d'un problème d'optimisation non-contraint sur des vecteurs de $\mathbb{R}^d$. Etant donné des réels $b_1,\ldots,b_n$ et des vecteurs $a_1,\ldots,a_n \in \mathbb{R}^d$, on considère le problème suivant:
$$\min_{x \in \mathbb{R}^d} \underbrace{-\sum^n_{i=1} \log(b_i-a_i^Tx)}_{≝ \; f(x)}$$
La résolution de ce problème s'appelle la recherche du centre analytique défini par les inégalités $a_i^Tx\leq b_i$
On pose
$$f ≝ \mathbb{R}^d \ni x ⟼ -\sum^n_{i=1} \log(b_i- ⟨a_i, x⟩)$$
NB: Si on suppose que les $a_i$ et $x$ sont des vecteurs lignes (ce qu'on fera dans la suite), $⟨a_i, x⟩$ vaut plutôt $a_i x^T$ (c'est un scalaire)
On a:
$$\begin{cases} \nabla f (x) &≝ \, \displaystyle \sum^n_{i=1} \frac{a_i}{b_i - ⟨a_i, x⟩} \\ \nabla^2 f (x) &≝ \, \displaystyle \sum^n_{i=1} \frac{\overbrace{a_i^T a_i}^{\text{matrice} \, (n,n)}}{(b_i-⟨a_i, x⟩)^2} \\ \end{cases}$$
qu'on peut réécrire sous la forme
$$\begin{cases} \nabla f (x) &≝ \; (A^T c(x))^T \\ \nabla^2 f (x) &≝ \; A^T C(x)^2 A \end{cases}$$
où
$c(x)$ est le vecteur colonne $\begin{pmatrix} \frac{1}{b_1 - ⟨a_1, x⟩} \\ \vdots \\ \frac{1}{b_n - ⟨a_n, x⟩} \end{pmatrix}$
$C(x) ≝ \operatorname{diag} (c(x))$ est la matrice diagonale de diagonale $c(x)$: $$C(x) ≝ \begin{pmatrix} \frac{1}{b_1 - ⟨a_1, x⟩}&0&\ldots &0\\ 0& \frac{1}{b_2 - ⟨a_2, x⟩} &\ddots &\vdots \\ \vdots &\ddots &\ddots &0\\ 0&\ldots &0&\frac{1}{b_n - ⟨a_n, x⟩}\end{pmatrix}$$
L'expression de $\nabla^2 f (x)$ vient du fait que si $X ≝ (X_1 | ⋯ | X_d)$ est une matrice dont les colonnes sont notées $X_i$, et $L ≝ \begin{pmatrix}L_1 \\ \hline \vdots \\ \hline L_n \end{pmatrix}$ une matrice dont les lignes sont notées $L_i$:
$$XL = \sum\limits_{i =1 }^d X_i L_i$$
Il vient alors que pour tout $x ∈ \operatorname{dom} f$, $\nabla^2 f (x)$ est bien symétrique (clair), et définie positive dès que $\operatorname{rang}(A) = d ≤ n$, puisque $∀X∈ 𝔐_{d, 1}(ℝ)$:
$$\begin{align*} & \nabla^2 f (x) X = 0 \\ ⟺ \; & A^T C(x)^2 A X = 0\\ ⟹ \; & X^T A^T C(x)^2 A X = 0\\ ⟹ \; & (C(x)AX)^T C(x) A X = 0\\ ⟹ \; & \Vert C(x) A X \Vert_2^2 = 0\\ ⟹ \; & C(x) A X = 0\\ ⟹ \; & A X = 0 &&\text{( } C(x) \text{ inversible)}\\ ⟹ \; & X = 0 &&\text{(puisque } \operatorname{rang}(A) = \dim 𝔐_{d, 1}(ℝ) \text{)} \end{align*}$$
d'où la stricte convexité de $f$.
# /!\: j'utilise PYTHON3
from IPython.display import display, Markdown
# Génération des données
d = 20 # dimension
n = 100 # nombre de points
b_inf, b_sup = 0, 1 # b_inf <= b_i <= b_sup
assert b_inf >= 0 # pour le que le x initial soit dans dom(f) plus tard
a_inf, a_sup = -1, 1 # a_inf <= a_i <= a_sup
A = (a_sup - a_inf) * np.random.rand(n, d) + a_inf # matrice dont les lignes sont les a_i
b = (b_sup - b_inf) * np.random.rand(n) + b_inf # vecteur dont les coefficients sont les b_i
# Initialisation
eps = 1e-8 # critère d'arrêt
alpha, beta = 0.01, 0.5 # paramètres du backtracking line search
n_iter_max = 600
def newton(A=A, b=b, eps=eps, n_max=n_iter_max, alpha=alpha, beta=beta, f_star=None, values_f=False):
def f(x):
return -sum(np.log(b[i]- np.dot(A[i], x)) for i in range(n))
# erreur commise et valeurs de f
err = []
f_values = []
# le x initial doit être dans le domaine de f, i.e. A x < b
# Ici, prenons x=0, comme b est positif
x = np.zeros(d)
for _ in range(n_iter_max):
# c(x) = [1/(b_1 - a_1^T x), ..., 1/(b_n - a_n^T x)]
# de sorte que grad f (x) = A^T c(x)
# et Hess(f)(x) = A^T (C(x))^2 A, où C(x) est la matrice diagonale tq C[i, i] = c(x)[i]
c = np.array([1/(b[i]- np.dot(A[i], x)) for i in range(n)])
grad_f = np.dot(A.T, c)
hess_f = np.dot(np.dot(A.T, np.diag(c)**2), A)
# pas de Newton
Delta_x = -np.dot(np.linalg.inv(hess_f), grad_f)
# -lambda^2 = (grad f)^T Delta_x
# (avec les notations de http://www.di.ens.fr/~aspremon/PDF/ENS/InteriorPointENS.pdf page 20)
grad_Delta_x = np.dot(grad_f, Delta_x)
if abs(grad_Delta_x)/2 <= eps:
break
### Début Backtracking line search
t = 1
# tant que x + t Delta_x n'est pas dans le dom f,
# i.e. tant que A (x + t Delta_x) < b n'est pas vrai
while (np.dot(A, x + t*Delta_x) >= b).any():
t *= beta
# f(x + t Delta_x) (stocké dans err et/ou f_values à la fin)
f_x_tDelta_x = f(x + t*Delta_x)
# tant que f(x + t Delta_x) > f(x) + alpha t (grad f)^T Delta_x
while f_x_tDelta_x > f(x) + alpha*t*grad_Delta_x:
t *= beta
f_x_tDelta_x = f(x + t*Delta_x)
### Fin Backtracking line search
# mise à jour de la distance à l'optimal et des valeurs de f
if f_star is not None:
err.append(abs(f_x_tDelta_x - f_star))
if values_f:
f_values.append(f_x_tDelta_x)
x += t*Delta_x
if f_star is not None:
return (x, f_values, err) if values_f else (x, err)
else:
return (x, f_values) if values_f else x
x_star, f_values = newton(values_f=True)
display(Markdown("- $A \; ≝ \; $ {} \n\n\n"\
"- $b \; ≝ \; $ {} \n\n\n"\
"- $x^\star \; ≝ \; $ {} \n"
.format(A, b, x_star)))
# affichage des valeurs prises par f pendant l'exécution de l'algorithme
data = go.Scatter(x = list(range(1, len(f_values)+1)), y = f_values, mode = 'lines', name = 'valeurs de f')
data_plot = [data]
layout= go.Layout(
title= 'Méthode de Newton avec backtracking line search',
hovermode= 'closest',
xaxis= dict(
title= 'Itérations',
ticklen= 5,
zeroline= False,
gridwidth= 2,
),
yaxis=dict(
title= 'f',
ticklen= 5,
gridwidth= 2,
),
showlegend= True
)
plotly.offline.iplot(go.Figure(data=data_plot, layout=layout))
L'optimisation de ce problème est une première étape vers des méthodes très largement utilisées à l'heure actuelle dans le domaine de l'optimisation non contrainte : les méthodes de point intérieur. Elles sont traitées in extenso dans le livre de Boyd et Vandenberghe servant de référence à cette partie du cours.
def f(x):
return -sum(np.log(b[i]- np.dot(A[i], x)) for i in range(n))
# calcul de la valeur de l'objectif x_star au bout d'un très grand nombre d'itération
x_star = newton(eps = 1e-10, n_max=10000)
# erreur commise
_, err = newton(eps = 1e-10, n_max=10000, f_star = f(x_star))
# affichage de l'évolution de l'erreur
data = go.Scatter(x = list(range(1, len(err)+1)), y = err, mode = 'lines', name = '|f_t - f^*|')
data_plot = [data]
layout= go.Layout(
title= 'Méthode de Newton (échelle logarithmique)',
hovermode= 'closest',
xaxis= dict(
title= 'Itérations',
ticklen= 5,
zeroline= False,
gridwidth= 2,
),
yaxis=dict(
title= 'Erreur commise',
type='log',
ticklen= 5,
gridwidth= 2,
),
showlegend= True
)
plotly.offline.iplot(go.Figure(data=data_plot, layout=layout))