In [7]:
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>'
))

TP 4 : Optimisation Convexe

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

Exercice 2: Centre analytique d'inégalités par méthode de Newton

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$

1) Implémenter la résolution de ce problème par méthode de Newton. Vous penserez à ajuster le pas par backtracking line search.

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

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

In [39]:
# /!\: 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))
  • $A \; ≝ \; $ [[ 0.90403743 -0.7160882 -0.38858907 ..., -0.39432233 0.45093446 -0.86071238] [-0.4213948 -0.52281251 0.41244602 ..., -0.68046322 0.53264587 0.26277129] [ 0.31679542 0.65268307 -0.7944629 ..., -0.38053765 0.77466676 -0.19098143] ..., [-0.86882736 -0.07585041 -0.79463656 ..., -0.18710827 0.0013474 -0.51631114] [-0.58736212 -0.53654925 -0.39825347 ..., 0.49291561 0.85920366 -0.86290787] [-0.22924026 -0.04349808 -0.59841209 ..., -0.19687436 -0.70417651 -0.69272616]]
  • $b \; ≝ \; $ [ 0.23586216 0.71365924 0.55441179 0.2213286 0.69645818 0.24115019 0.13061103 0.92057883 0.83688339 0.46065425 0.33172266 0.08828906 0.27622202 0.1897234 0.65281579 0.12060597 0.62537309 0.18569229 0.06683669 0.06177858 0.45217718 0.60696985 0.22428507 0.58087361 0.16198375 0.93516555 0.62414533 0.17952969 0.08311998 0.57964939 0.01420207 0.58045763 0.0432707 0.43054385 0.16753945 0.03567748 0.77229234 0.4120146 0.97780193 0.30405315 0.89129666 0.21088253 0.32357774 0.46091643 0.02199559 0.35620951 0.41212106 0.18222654 0.19615615 0.75188682 0.82980847 0.59474721 0.35599823 0.07056911 0.01224795 0.98038105 0.21849597 0.24379085 0.55574603 0.57574261 0.75705027 0.81719679 0.67411482 0.65106481 0.40183226 0.48204334 0.65324783 0.58844588 0.81192283 0.51371338 0.90260994 0.79368228 0.94344192 0.81738762 0.03484231 0.51391057 0.07589396 0.49231416 0.59519131 0.67026785 0.37261569 0.94729369 0.09828982 0.22127713 0.07425842 0.12884968 0.78131552 0.64316873 0.50273615 0.69565892 0.09100803 0.52737543 0.71723139 0.58674221 0.54995826 0.23471469 0.73216917 0.97992772 0.05031723 0.32377069]
  • $x^\star \; ≝ \; $ [ 0.08748678 -0.08918145 0.01651846 -0.05297649 0.08288545 -0.10638417 0.00503012 0.02499581 -0.16047685 0.04710907 0.05015083 0.01197512 0.01856327 -0.09142799 0.17961693 -0.04922944 0.06000178 0.00041643 -0.15040295 -0.03072997]

2) Essayez de représenter graphiquement l'évolution de l'erreur commise $\|f_t-f^*\|_2$ en fonction du nombre d'itérations de l'algorithme(on pourra pour cela utiliser une méthode très simple: calculer la valeur de l'objectif au bout d'un très grand nombre d'itération et prendre cette valeur comme approximation de la vraie valeur de la fonction à l'optimum $f^*$).

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.

In [41]:
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))