Algorithmes d'optimisation -- L3 MINT et doubles licences 2019/2020 -- 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}}$
Dans ce TP, on cherche à appliquer l'algorithme d'Uzawa au calcul de la projection d'un point $p\in \Rsp^d$ sur un polyèdre $K$, c'est-à-dire un ensemble convexe défini par un nombre fini d'inégalités affines:
$$ 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$. Dans la suite, on note $x\leq y$ où $x,y$ sont deux vecteurs si $\forall i, x_i\leq y_i$. En posant $A$ la matrice possédant $k$ lignes notées $a_1,\hdots,a_k \in \Rsp^d$ et $b\in\Rsp^k$, on a donc
$$ K = \{ x \in \Rsp^d\mid Ax \leq b\} $$Dans cette première partie, on donne un algorithme permettant de calculer la projection d'un point $p \in \Rsp^d$ sur le polyèdre $K$:
$$ (P) := \min_{x\in K} \frac{1}{2} \nr{x - p}^2 $$Le lagrangien $L$ du problème (P) est donné par
$$L: (x,\lambda)\in\Rsp^d\times \Rsp^k_+ \mapsto f(x) + \sum_{1\leq i\leq k} \lambda_i c_i(x)$$où $f(x) = \frac{1}{2} \nr{x - p}^2$, et le problème dual associé à (P) est donc
$$ (D) := \max_{\lambda \in \Rsp^k_+} \min_{x\in \Rsp^d} L(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} L(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.
En particulier, le problème dual est un problème d'optimisation avec contraintes ($x\in M$), mais l'ensemble de contraintes est très simple.
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)})$.
Ainsi, si la suite $(\lambda^{(k)})$ converge, sa limite est un maximiseur $\lambda^*$ de (D) et $x_{\lambda^*} = x^*$ est un minimiseur de (P).
Q3) [Mise en oeuvre] Mettre en oeuvre l'algorithme d'Uzawa.
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)}$.
projection_convexe(A,b,p,tau)
, pour un assez grand nombre (100) de points p choisis aléatoirement dans $[-4,4]^2$.import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# la commande suivante agrandit les figures
plt.rcParams['figure.figsize'] = [9.,6.]
# <completer>
# <completer>
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
À 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.
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.
p = t**4 + .2*np.random.rand(n)
# <completer>