Algorithmes d'optimisation -- L3 MINT et doubles licences 2017/2018 -- Université Paris-Sud


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

TP 5: Projection, débruitage et algorithme d'Uzawa

 Partie I: Étude et mise en oeuvre de l'algorithme d'Uzawa

Dans ce TP, on cherche à appliquer l'algorithme d'Uzawa au calcul de la projection d'un point $p\in \Rsp^d$ sur un convexe $K$ défini par des inégalités affines. Dans la suite, on note $x\leq y$ où $x,y$ sont deux vecteurs si $\forall i, x_i\leq y_i$.

  • On considère $K = \{ x \in \Rsp^d\mid Ax \leq b\}$, où $A$ est une matrice possédant $k$ lignes notées $a_1,\hdots,a_k \in \Rsp^d$ et $b\in\Rsp^k$. De manière équivalente, $$ K = \{ x\in \Rsp^d \mid \forall 1\leq i\leq k, c_i(x) \leq 0 \} $$ où $c_i(x) := \sca{a_i}{x} - b_i \leq 0$. On s'intéresse au problème de projection d'un point $p \in \Rsp^d$ sur $K$: $$ (P) := \min_{x\in K} \frac{1}{2} \nr{x - p}^2 $$
  • Le lagrangien $\ell$ du problème (P) est donné par $\ell: (x,\lambda)\in\Rsp^d\times \Rsp^k_+ \mapsto f(x) + \sum_{1\leq i\leq k} \lambda_i c_i(x)$ et le problème dual associé à (P) est donc $$ (D) := \max_{\lambda \in \Rsp^k_+} \min_{x\in \Rsp^d} \ell(x,\lambda) $$
  • On rappelle que si $\lambda^*$ est un maximiseur de (D), alors tout solution du problème de minimisation sans contrainte $(P_{\lambda^*}) = \min_{x\in\Rsp^d} \ell(x,\lambda^*)$ est aussi solution du problème (P). En d'autre terme, la connaissance de $\lambda^*$ permet de remplacer un problème d'optimisation avec contraintes $(P)$ par un problème d'optimisation sans contrainte $(P_{\lambda^*})$ !

Nous allons étudier dans cette partie l'algorithme d'Uzawa. L'idée est de calculer un maximiseur $\lambda^*$ du problème dual (D) par une méthode de gradient projeté, et de s'en servir pour calculer la solution de (P) en utilisant le dernier rappel.

Q1) [Expression du problème dual] Dans cette question, il s'agit d'écrire le problème dual de manière plus explicite.

  • Montrer que le lagrangien associé au problème (P) peut s'écrire $\ell(x,\lambda) = \frac{1}{2}\nr{x - p}^2 + \sca{\lambda}{A x - b}$.
  • Étant donné $\lambda \in \Rsp^k$, donner l'expression de l'unique solution $x_\lambda$ du problème de minimisation $\min_{x \in \Rsp^d} \ell(x,\lambda).$
  • En déduire que $$\begin{aligned} &\qquad (D) := \max_{\lambda \in\Rsp_+^k} -h(\lambda) \\ &\hbox{ où } h(\lambda) = \frac{1}{2} \nr{A^T \lambda - p}^2 - \frac{1}{2}\nr{p}^2 + \sca{\lambda}{b} \end{aligned}$$
  • Montrer que $\nabla h(\lambda) = A(A^T \lambda - p) + b = b - A x_\lambda$.

Algorithme d'Uzawa: on appelle ainsi l'algorithme du gradient projeté pour le problème dual (D): $$ \begin{cases} \lambda^{(0)} = 0 \in \Rsp^k \\ g^{(k)} = \nabla h(\lambda^{(k)})\\ \lambda^{(k+1)} = p_{\Rsp_+^k}(\lambda^{(k)} - \tau g^{(k)})\\ x^{(k+1)} = p - A^T \lambda^{(k+1)} \quad (\in \arg\min_{x\in\Rsp^d} \ell(x,\lambda^{(k+1)})) \end{cases} $$ L'algorithme est arrêté lorsque $\nr{x^{(k)} - x^{(k+1)}}\leq \eps$.

Pour l'implémentation de l'algorithme, on rappelle que $p_{\Rsp_+^k}(v) = (\max(v_1,0),\hdots,\max(v_k,0))$.

Q2) [Convergence de l'algorithme d'Uzawa] On pose $S_\tau(\lambda) := p_{\Rsp_+^k}(\lambda - \tau \nabla h(\lambda))$, de sorte que $\lambda^{(k+1)} = S_\tau(\lambda^{(k)})$.

  • Montrer que la fonction $h$ est convexe.
  • En déduire que si $\lambda^*$ est un point fixe de $S_\tau$, alors $\lambda^*$ est solution du problème (D) (i.e. maximise $h$ sur $\Rsp_+^k$).

Ainsi, si la suite $(\lambda^{(k)})$ converge, sa limite est un maximiseur $\lambda^*$ de (D) et $x_{\lambda^*} = x^*$ est un minimiseur de (P).

Q4) [Mise en oeuvre] Mettre en oeuvre l'algorithme d'Uzawa.

  • Écrire une une fonction projection_convexe(A,b,p,tau,err=1e-6) calculant les itérées de $(\lambda^{(k)}, x^{(k)})$, en arrêtant la boucle dès que $\nr{x^{(k)}- x^{(k+1)}} \leq$ err. Cette fonction retournera $\lambda^{(k)}, x^{(k)}$.
  • Tester cette fonction sur le convexe $K = \{ x \in \Rsp^2\mid \nr{x}_1 \leq 1 \}$. On commencera par déterminer $A,b$ tel que $K = \{x \mid Ax \leq b \}$. On vérifiera la validité du calcul de deux manières:
    • visuellement, en affichant le segment reliant p à son projeté q = projection_convexe(A,b,p,tau), pour un assez grand nombre (100) de points p choisis aléatoirement dans $[-4,4]^2$.
    • en vérifiant que la solution $x,\lambda$ = projection_convexe(A,b,p,tau) satisfait (à erreur numérique près) les quatre conditions du théorème de Karush-Kuhn-Tucker: $Ax \leq b$ (admissibilité de $x$), $\lambda \geq 0$ (admissibilité de $\lambda$), $\forall i, (A x - b)_i \lambda_i = 0$ (complémentarité) et $\nabla_x \ell(x,\lambda) = 0$ (optimalité).
  • Recommence avec $K = \{ x \in\Rsp^2 \mid \nr{x}_\infty \leq 1 \}$.
In [5]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# la commande suivante agrandit les figures
plt.rcParams['figure.figsize'] = [9.,6.]

# <completer>
In [6]:
# <completer>

Partie II: Régression isotone

Nous allons considérer deux problèmes de débruitage consistant simplement à projeter sur un convexe. Les données sont par exemple des séries temporelles $y = (y_1,\hdots,y_n)\in \Rsp^n$, mesurées avec un bruit. On sait que les données réelles appartiennent à un certain ensemble convexe $K$ de $\Rsp^n$. Deux exemples

  • régression isotone: $K = \{ x\in \Rsp^n \mid \forall 1\leq i < n,~x_{i+1}\geq x_i \}$
  • régression convexe: $K = \{ x\in \Rsp^n \mid \forall 1 < i < n,~x_{i} \leq \frac{1}{2} (x_{i-1} + x_{i+1}) \}$.

À cause du bruit, le vecteur $y$ mesuré n'appartient pas à l'ensemble $K$. L'idée est simplement de débruiter le signal en le projetant sur $K$, soit:

$$ (P) := \min_{x\in K} \frac{1}{2} \nr{x - p}^2 $$

Q1) Implémenter la régression isotone.

  • Trouver une matrice $A$ et un vecteur $b$ tel que $K_{iso} = \{ x\in \Rsp^n \mid Ax \leq b \}$.
  • Utiliser projection_convexe() avec tau=0.1 et avec le vecteur $p$ donné ci-dessous.
  • Que peut-on dire expérimentalement de la solution $x^*$ aux points $\{i,i+1\}$ où $y_{i}\geq y_{i+1}$ ?
  • Démontrer ce résultat en utilisant les conditions du théorème KKT.
In [8]:
n = 30
t = np.linspace(0,1,n)
p = t**2 + .3*np.random.rand(n)

# <completer>

Q2) Implémenter la régression convexe en utilisant projection_convexe() avec tau=0.1 et avec le vecteur $p$ donné ci-dessous.

In [9]:
p = t**4 + .2*np.random.rand(n)

# <completer>