\chapter{Optimisation pour l'analyse de données}
\label{ch:optimisation}

%======================================================================
% MY MACROS
%======================================================================

% dans bookDS
%\DeclareMathOperator{\prox}{prox}

\def\RR{\mathbb{R}}
\def\eps{\varepsilon}
\def\si{\sigma}
\def\al{\alpha}
\def\la{\lambda}
\def\ga{\gamma}
\def\Ph{\Phi}
\def\w{{\pmb \omega}}
\def\y{\mathbf{y}}
\def\z{\mathbf{z}}
\def\u{\mathbf{u}}
\def\x{\mathbf{x}}
\def\v{\mathbf{v}}
\def\a{\mathbf{a}}
\def\om{\mathbf{\omega}}
%\def\om{\mathbf{\omega}}
\def\Om{\Omega}
%\DeclareMathOperator{\prox}{prox}
\newcommand{\ip}[2]{\langle #1, #2 \rangle}
%\newcommand{\norm}[1]{\left\| #1 \right\|}
\newcommand{\cad}{c'est-à-dire\xspace}
\newcommand{\trans}[1]{{#1}^{\top}}
\newcommand{\nom}[1]{{\emph{#1}}}
\newcommand{\argmin}{\text{argmin}}
%\newcommand{\guill}[1]{{"{#1}}"}


\newcommand{\accol}[1]{{\left\{ \begin{array}{ll} #1 \end{array} \right.}}


\newcommand{\illustrationtextscale}[4]{
    \noindent
    \begin{minipage}[c]{#1\linewidth}
      #3
    \end{minipage} \hfill
    \begin{minipage}[c]{#2\linewidth}
      \includegraphics[width=\linewidth]{#4}
    \end{minipage}
}




%======================================================================
%======================================================================



L'optimisation est une discipline des mathématiques appliquées qui 
se retrouve au c{\oe}ur des méthodes numériques pour l'analyse de données. Par exemple, en apprentissage statistique, des algorithmes d'optimisation exploitent de grands volumes de données pour calculer les paramètres des modèles utilisés ensuite pour prédire, classifier ou décider. 

%======================================================================
\paragraph{Objectifs de ce chapitre}

Ce chapitre présente une vue d'ensemble de l'optimisation mathématique dans le contexte des applications en science des données. L'optimisation a plusieurs aspects: la formulation des problèmes eux-mêmes, la conception d'algorithmes pour les résoudre, l'étude de leurs propriétés théoriques et leur mise en {\oe}uvre pratique. Chacun de ces points fait l'objet de livres entiers et de nombreux articles de recherche; ici, nous n'allons qu'effleurer ces questions, dans la perspective de l'analyse de données. Nous allons insister sur la structure particulière des problèmes d'optimisation apparaissant en recommandation, classification et régression, et impliquant des données de grande taille, éventuellement stockées sur différentes machines. Nous rappellerons les méthodes de base de l'optimisation et les idées qui permettent aux algorithmes de passer à l'échelle. Ce chapitre est écrit à destination d'étudiants suivant un master\footnote{Les pré-requis sont uniquement des notions de niveau bac+3 en mathématiques et informatique. Nous manipulons des objets basiques en mathématiques (matrices, fonctions, gradients), les notions basiques des probabilités et statistiques (espérance, variance, moyenne empirique) et des concepts algorithmiques simples (itération, coût de calcul, critère d'arrêt).}
 à coloration \guil{science des données} et aux étudiants en thèse souhaitant une introduction à l'optimisation dans ce cadre particulier.  

%======================================================================
\paragraph{Trois aspects de l'optimisation en science des données}

Comment l'optimisation appara\^{\i}t-elle en analyse des données ? Quelles sont les méthodes pour résoudre les problèmes d'optimisation de grande taille ? Ce chapitre apporte ainsi des éléments de réponse, en insistant sur trois aspects complémentaires, particulièrement intéressants dans un contexte de méga-données:
\begin{itemize}
\item les algorithmes du premier ordre qui réalisent un bon compromis entre simplicité et efficacité;%permettent d'obtenir rapidement des solutions numériques de précision moyenne;
\item les algorithmes aléatoires qui traitent des sous-parties des données;
\item le calcul parallèle et distribué, pour lequel les méthodes de premier ordre sont bien adaptées.
\end{itemize}
%Nous illustrerons ces aspects avec des problèmes classiques apparaissant en recommendation, classification et régression.

%======================================================================
%======================================================================
\section{Apprentissage et optimisation}\label{sec:gene}
%======================================================================
%======================================================================
%L'optimisation est un pilier de l'apprentissage automatique. Le calcul des paramètres d'un système qui prend des décisions sur de nouvelles données. Passer sur les données disponibles ces paramètres sont choisis \nom{au mieux} pour le modèle d'apprentissage choisi. 

La cha\^{\i}ne de traitement des données se résume en trois grandes étapes:
observer le monde (collecter des données),
proposer des modèles (conception et apprentissage),
tester sur de nouvelles données (estimer l'erreur). %, appelée erreur de généralisation)
La phase d'apprentissage consiste à chercher les paramètres d'un modèle stochastique expliquant \nom{au mieux} les données et c'est ici qu'apparaît l'optimisation, lors du calcul des paramètres optimaux. Nous illustrons ceci dans le contexte général de l'apprentissage supervisé; l'apprentissage non supervisé\index{apprentissage!non-supervisé} ou semi-supervisé\index{apprentissage!semi-supervisé} (ou même d'autres techniques d'apprentissage) amène de la même manière à des problèmes d'optimisation.



%======================================================================
\paragraph{Notations, définitions}

Commençons par définir quelques termes et notations.  Dans ce chapitre, la variable d'optimisation est un vecteur de taille $d$ qui sera noté $\w\in \RR^d$, en cohérence avec le reste du livre\footnote{Dans les livres et les articles d'optimisation, il est coutume de noter la variable\;$\x$; en statistique, c'est plutôt $\beta$. Dans ce livre et ce chapitre en particulier, nous suivons la notation standard en apprentissage: la variable du problème d'optimisation, correspondant au paramètre du modèle d'apprentissage, est souvent notée $\w\in \RR^d$.}.
On note le produit scalaire canonique de l'espace $\RR^d$ (somme des produits terme à terme de deux vecteurs):
\[
\ip{\w}{\v}=\sum^d_{j=1} \om_jv_j \qquad\text{pour $\w\in \RR^d$ et $\v\in \RR^d$}
\]
Dans ce chapitre, on utilise deux normes dans $\RR^d$: la norme\;$\ell_2$ (la norme euclidienne associée au produit scalaire) et la norme\;$\ell_1$, définies respectivement par :
\[
\norm{\w}_2=\sqrt{\ip{\w}{\w}}%=\sqrt{\sum^d_{i=1} {x_i}^2}
=\Big(\sum^d_{j=1} {\om_j}^2\Big)^{\frac{1}{2}}
\qquad\text{et}\qquad
\norm{\w}_1=\sum^d_{j=1} |\om_j | \qquad\text{pour $\w\in \RR^d$}
\]

Pour les problèmes d'optimisation qui font plus particulièrement l'objet de la section \ref{sec:optim}, nous utilisons la terminologie suivante.

On considère une fonction $f\colon\RR^d\rightarrow\RR$, qu'on appelle \nom{fonction-objectif} ou \nom{critère}, et un ensemble $C\subset \RR^d$, qu'on appelle \nom{ensemble des contraintes} ou \nom{ensemble réalisable}. Le problème est alors de trouver la valeur optimale $f_\star\in \RR$ et une solution optimale $\w_{\star}\in C$ telles que:\index{solution optimale}\index{contrainte}\index{objectif}
\[
f_{\star} =f(\w_{\star})~ \leq ~f(\w) \qquad\text{pour tout } \w\in C\,
\]
Ceci s'écrit formellement sous la forme:\index{problème d'optimisation}
\begin{equation}\label{eq:opt}
\min_{\w\in C} ~~f(\w)
\end{equation}


%en considérant des problèmes linéaires et des problèmes non linéaire et non convexe en apprentissage profond.
%======================================================================
\paragraph{Contexte d'apprentissage supervisé}

Pla\c{c}ons-nous dans un cadre standard d'apprentissage supervisé.\index{apprentissage!supervisé}\index{supervised learning|see{apprentissage supervisé}}
Après une phase préliminaire de collecte, stockage et nettoyage des données, nous disposons d'un ensemble de données d'entra\^{\i}nement constitué de $n$ observations sous forme de couples $(\x_i,y_i)$ avec un vecteur de \nom{caractéristique} $\x_i\in \RR^m$ et un \nom{attribut}  $y_i\in \RR$ (appelé aussi \nom{étiquette}, \nom{score} ou \nom{label}). Pour fixer les idées, on peut penser que $\x_i$ représente les réponses d'un étudiant à un questionnaire sur sa vie personnelle et $y_i$ représente sa note à l'examen final (ainsi $y_i\in [0,20]$) ou alors l'information qu'il valide ou non son année (dans ce cas, $y_i$ est binaire). 
%Dans cette section, on note ainsi $n$ taille de l'échantillon et $D$ la taille des caractéristiques.
Nous souhaitons alors utiliser ces données pour \emph{prédire} l'étiquette associée à une nouvelle donnée de caractéristique $a$. Dans l'exemple des notes des étudiants, on voudrait avoir, à partir de celles des étudiants de l'année précédente, une prédiction des notes finales des nouveaux étudiants. En général, le problème de prédiction s'appelle \nom{régression}  dans le cas où $y_i\in \RR$ et \nom{classification} dans le cas où les $y_i$ prennent un nombre fini de valeurs (par exemple, classification binaire si $y_i\in \{-1,1\}$).\index{régression}\index{classification}\index{prediction@prédiction}

En apprentissage supervisé, cette prédiction se fait par un modèle statistique calibré sur les données d'entra{\^{\i}}nement.
On choisit généralement un modèle avec une forme fixée dépendante d'un vecteur de paramètres que l'on note $\w\in\RR^d$. Une fonction de prédiction générale est ainsi une fonction $h(\cdot,\w)$ paramétrée par $\w$ qui propose une prédiction $b=h(\x,\w)\in \RR$ pour une nouvelle donnée $\x$.
Voici deux exemples fondamentaux de fonctions de prédiction:
\begin{itemize}
\item Prédiction linéaire : $h(\x,\w) = \ip{\x}{\w}$ dans le cas le plus simple (ou $h(\x,\w) = \ip{\phi(\x)}{\w}$ dans le cas d'un espace d'attributs transformés pour introduire de la non-linéarité); 
\item Prédiction (hautement) non linéaire par réseaux de neurones artificiels : $h(\x,\w) = \ip{\w_m}{\si(\ip{\w_{m-1}}{\cdots\si(\ip{\w_1}{\x})})}$ 
où la prédiction est calculée en appliquant une succession de compositions de combinaisons linéaires et d'une fonction d'activation non linéaire (comme la fonction sigmo\"{\i}de $\si(z)=1/(1+\exp(-z))$ ou la fonction ReLU $\si(z)=\max\{0,z\}$).
\end{itemize}

%– If fi(θ) = l(yi, Φ(xi)⊤θ), then fi′(θ) = l′(yi, Φ(xi)⊤θ) Φ(xi)

%======================================================================
\paragraph{Apprendre, c'est optimiser}

La phase d'apprentissage consiste à calculer le paramètre $\w\in \RR^d$ du modèle 
qui explique \nom{au mieux} les données disponibles. Cela signifie que la prédiction doit 
\nom{être proche} des vraies étiquettes sur les données connues :
\begin{equation}\label{eq:train}
h(\x_i,\w) \simeq y_i\qquad \text{pour $(\x_i,y_i)$ données d'entra{\^{\i}}nement}
\end{equation}
et aussi pouvoir bien se généraliser pour des données inconnues.
On mesure la différence \eqref{eq:train} entre la prédiction et l'étiquette avec une fonction, appelée \nom{fonction perte}, 
$\ell\colon\RR\times\RR\rightarrow\RR$, dont la plus simple et connue est la perte des moindres carrées $\ell(y,z) = \frac{1}{2}(y-z)^2$. D'autres fonctions pertes sont très utilisées aussi, comme $\ell(y,z) = |y-z|$ en régression robuste, 
$\ell(y,z) = \log(1+ \exp(-y\,z))$ en régression logistique, ou
$\ell(y,z) = \max\{0,1-y\,z)$ en classification.
Il est naturel de s'intéresser alors au paramètre qui réalise le moins d'erreurs \eqref{eq:train} au sens de la mesure choisie $\ell$, \cad $\w\in \RR^d$ solution du problème d'optimisation:
\begin{equation}\label{eq:minERM0}
\min_{\w\in \RR^d} ~~ \frac{1}{n}\sum^n_{i=1} \ell\big(y_i, h(\x_i,\w)\big)
\end{equation}


%======================================================================
\paragraph{Point de vue statistique et conséquences pratiques} 

Le problème \eqref{eq:minERM0} admet une interprétation statistique dans le cas où les données sont des réalisations d'une variable aléatoire. En effet, si on suppose que les données $(\x_i,y_i)$ sont un échantillon des observations indépendantes et identiquement distribuées d'un vecteur aléatoire $(\x,y)\in \RR^{d+1}$, le problème d'optimisation sous-jacent consiste en la minimisation de l'espérance de l'erreur, qui s'écrit :
\begin{equation}\label{eq:minRM}
\min_{\w\in \RR^d} ~~\EE\left[\ell\big(y,h(\x,\w)\big)\right]
\end{equation}
où l'espérance est prise sur la distribution de $(\x,y)$. % la mesure de probabilité $\PP$ de $(\x,y)$. 
La perte peut aussi être une estimation de maximum de vraisemblance %$-\log \PP_\w\big((\x,y_i)\big)$ 
pour une distribution paramétrée par $\w$.
La fonction-objectif \eqref{eq:minERM0} s'interprète alors simplement comme l'estimateur empirique de l'espérance de l'erreur \eqref{eq:minRM}. En théorie, par la loi des grands nombres, l'espérance empirique va converger en loi vers l'espérance lorsque le nombre de données grandit. En pratique néanmoins, les données sont en nombre fini et l'hypothèse d'indépendance sous-jacente est discutable: nous travaillons donc directement sur l'espérance empirique \eqref{eq:minERM0} plutôt que de considérer \eqref{eq:minRM}. 

Il faut retenir de cette interprétation statistique que, dans une situation idéale, la minimisation de \eqref{eq:minERM0} n'est pas une fin en soit car elle néglige l'erreur statistique (la différence entre les valeurs optimales de \eqref{eq:minERM0} et\;\eqref{eq:minRM}). Il faut en fait avoir des mécanismes de contrôle de \nom{l'erreur totale}, qui est la somme des deux erreurs: l'erreur d'optimisation (la différence entre la solution optimale \eqref{eq:minERM0} et l'itéré courant d'un algorithme d'optimisation) et l'erreur statistique. Cette remarque incite tout d'abord à ne pas minimiser l'erreur d'optimisation trop fortement au risque de détériorer l'erreur statistique; cela justifie alors de choisir des algorithmes, comme ceux de la section\;\ref{sec:algo}, obtenant rapidement des solutions de précision moyenne. Cette remarque amène aussi à l'introduction d'un terme de régularisation dans le problème d'optimisation, comme présenté ci-après. 


%======================================================================
\paragraph{Minimisation du risque empirique régularisé} 


Plutôt que de considérer directement \eqref{eq:minERM0}, il est préférable d'ajouter dans l'objectif une fonction de régularisation  $r\colon\RR^d\rightarrow\RR$ pondérée par un paramètre $\la>0$, ce qui amène au problème d'optimisation:\index{régularisation}
\begin{equation}\label{eq:minERM}
\min_{\w\in \RR^d}~~ f(x) ~=~ \frac{1}{n}\sum^n_{i=1} \ell\big(y_i, h(\x_i,\w)\big) ~+~ \la \, r(\w)
\end{equation}
Par exemple, dans le cas d'un modèle linéaire $h(\x,\w) = \ip{\x}{\w}$, d'une perte par moindres carrés et d'une régularisation par norme $2$ (appelée aussi régularisation de Tikhonov), l'apprentissage revient au problème :
\begin{equation}\label{eq:ridge}
\min_{\w\in \RR^d} ~~\frac{1}{n} \sum^n_{i=1}\frac{1}{2}(y_i-\ip{\x_i}{\w})^{2} + \frac{\la}{2}{\norm{\w}_2}^{\!2} 
~ = ~ \frac{1}{2n}{\norm{X\w-\y}_2}^{\!2} + \frac{\la}{2}{\norm{\w}_2}^{\!2}
\end{equation}
où on rassemble toutes les données dans la matrice $X\in \RR^{n\times d}$ dont les lignes sont les $\trans{\x_i}$ et le vecteur
$\y\in \RR^n$ des $y_i$. Insistons sur le fait que le paramètre de régularisation $\la>0$ est un paramètre du \emph{problème d'optimisation}, mais pas un paramètre du \emph{modèle d'apprentissage} comme $\w\in \RR^d$.  

On remarque les deux choix extrêmes: avec $\la$ proche de $0$, on souhaite coller aux données d'entra\^{\i}nement, alors qu'à l'inverse, $\la$ très grand permet de diminuer l'importance des données. Entre ces deux choix, la présence de la régularisation $r$ a plusieurs intérêts:
\begin{itemize}
\item apprendre de manière plus raisonnée : éviter le surapprentissage sur les données connues pour mieux généraliser aux nouvelles données;\index{surapprentissage}\index{overfitting|see{surapprentissage}}
\item promouvoir une structure particulière sur le paramètre d'apprentissage, par exemple de la parcimonie 
avec $r(\om)=\norm{\om}_1$;
\item améliorer la résolution numérique du problème d'optimisation, par exemple en améliorant le conditionnement avec la norme $r(\om)={\norm{\om}_2}^2$ (voir section \ref{sec:algo}). 
\end{itemize}

%======================================================================
\paragraph{Différentiabilité et gradient}

Terminons cette formalisation des problèmes d'optimisation en apprentissage par une remarque importante en vue de la résolution numérique. Certaines modélisations sous la forme \eqref{eq:minERM} considèrent un modèle $h$, une perte~$\ell$ et une régularisation $r$, tous trois différentiables, ce qui implique que la fonction objectif $f$ est différentiable aussi. Souvent,  le gradient $\nabla f(x)$ est de plus \emph{facilement calculable}, ce qui signifie que l'on dispose d'une expression explicite simple du gradient\index{gradient}, ou d'un algorithme efficace de calcul de dérivées itératives (appelée \guil{différentiation automatique} en optimisation ou \guil{propagation inverse} en apprentissage). Par exemple, la fonction objective du problème 
de régression ridge \eqref{eq:ridge} est différentiable et on a:
\begin{equation}\label{eq:gradridge}
\nabla f (\w)= \left(\frac{1}{n}\trans{X}\!X + \la \Id\right)\w - \frac{1}{n}\trans{X}\y
\end{equation}
%dont le coût de calcul est dominé par le coût de la multiplication de la matrice $X$ et sa transposée $\trans{X}$.
Cependant les fonctions objectives des problèmes d'apprentissage ne sont pas toujours différentiables: en fonction de la nature des données, la modélisation peut privilégier des pertes non différentiables ou des régularisations non différentiables; donnons deux exemples classiques: 
\begin{itemize}
\item la classification par \guil{SVM} (support vector machines) dont la perte admet une cassure à cause du max:
\begin{equation}\label{eq:SVM}
\min_{\w\in \RR^d} ~~ \frac{1}{n}\sum^n_{i=1} \max\{0,1-y_i\trans{\x_i}\w\}  ~+~ \frac{\la}{2}{\norm{\w}_2}^{\!2}
\end{equation}

\item la régression \guil{lasso} utilisant la norme $\norm{\cdot}_1$ pour induire des solutions optimales parcimonieuses (\cad avec beaucoup de $\w_i=0$), mais qui est justement non différentiable partout: \index{parcimonie} %quand une composant $\w_i$ s'annule
\begin{equation}\label{eq:lasso}
\min_{\w\in \RR^d} ~~ \frac{1}{2n}{\norm{X\w-\y}_2}^{\!2} ~+~ \la\norm{\w}_1\,
\end{equation}
\end{itemize}


%======================================================================
\paragraph{Bilan sur les problèmes d'optimisation en science des données}

Les problèmes d'optimisation en science des données sont typiquement de la forme :
\begin{equation}\label{eq:sumfinie}
\min_{\w\in \RR^d} ~~ \frac{1}{n}\sum^n_{i=1} f_i(\w) + r(\w)
\end{equation}
La minimisation se fait sur $\w$, le paramètre du modèle que l'on souhaite calculer; la fonction de régularisation $r:\RR^d\rightarrow\RR$ permet de \guil{mieux} apprendre; les fonctions $f_i:\RR^d\rightarrow\RR$ mesurent la qualité de prédiction associée à $\w$ sur un bloc de données. 
C'est ce que nous venons de voir, dans un cadre d'apprentissage supervisé, avec la minimisation du risque empirique régularisé \eqref{eq:minERM}. La résolution de ces problèmes d'optimisation forme la motivation et le c{\oe}ur de ce chapitre. 

En général, la difficulté particulière des problèmes d'optimisation issus de l'analyse de données ne vient pas de la formulation \eqref{eq:sumfinie} en elle-même, mais de la grande taille: grands volumes de données ($n$ grand, \emph{big data}) et/ou grande complexité des modèles ($d$ grand ou $h$ complexe, \emph{big models}). 
Ainsi la résolution de \eqref{eq:sumfinie} avec $n$ et $d$ fait l'objet des sections suivantes: nous allons tout d'abord revoir, dans la section \ref{sec:optim}, les idées de base de l'optimisation, puis étudier, dans la section\;\ref{sec:appli}, comment elles s'adaptent pour minimiser de grandes sommes finies comme \eqref{eq:sumfinie}. 

%Comment résoudre ce problème quand $d$ ou $n$ sont grands ? 



%======================================================================
%======================================================================
\section{Introduction à l'optimisation}\label{sec:optim}
%======================================================================
%======================================================================

Cette section donne une vue globale sur les fondamentaux de l'optimisation. La présentation insiste sur les idées en cachant sous le tapis les détails mathématiques et algorithmiques; l'objectif est uniquement de donner des éléments pour mieux comprendre les méthodes utilisées dans l'analyse des données. Un lecteur ayant déjà suivi un cours d'optimisation peut aisément sauter cette section.
%Utilisation de $x$ pour la prédiction, comme expliqué à la section précédente.
%On ne cherche pas de solutions avec une grande précision car les modèles sont inexacts.

%======================================================================
\subsection{Problèmes d'optimisation}
%======================================================================

Au-delà de son rôle en analyse de données, l'optimisation appara{\^\i}t dans de nombreuses disciplines scientifiques et dans les sciences de l'ingénieur au sens large (industrie, informatique, services).
Insistons cependant sur le fait que les problèmes apparaissent rarement déjà formalisés: une étape de modélisation est primordiale, notamment en analyse de données, comme illustré à la section précédente. Dans cette section, nous nous plaçons après cette étape de modélisation et nous considérons un problème d'optimisation, sous la forme générale \eqref{eq:opt}. Commençons par rappeler quelques idées à garder à l'esprit.


%======================================================================
\paragraph{Résoudre un problème d'optimisation}

Il est extrêmement rare que l'on puisse trouver une expression explicite des solutions optimales $f_{\star}$ et $\w_{\star}$. Dans quelques cas très simples, on peut tout de même en obtenir en résolvant explicitement des équations satisfaites par ces solutions. On les appelle les \emph{conditions d'optimalité}, qui sont, dans le cas particulier de $f$ différentiable sans contrainte:\index{condition!d'optimalité}
\begin{equation}\label{eq:condopt}
\text{$\w$ est un minimum (local) de $f$ différentiable ~$\Longrightarrow$~ $\nabla f(\w) =0$}
\end{equation}
En utilisant cette propriété, on montre par exemple que la solution optimale de \eqref{eq:ridge} est:
\[
\w_{\star} = -(\trans{X}X + \la \Id)^{-1}\trans{X} \y 
\]
et peut ainsi se calculer par la résolution du système linéaire.

Dans la plupart des cas, on approche donc numériquement $f_{\star}$ et $\w_{\star}$ par un algorithme itératif (voir section\;\ref{sec:algo}). On cherche ainsi une solution approchée à une précision $\eps$; cela signifie que l'on cherche $\bar \w\in C$ telle que $\norm{\bar \w - \w_{\star}} \leq \eps$, ou $f(\bar \w) - f_{\star} \leq \eps$, ou bien que les conditions d'optimalité sont vérifiées à $\eps$ près (par exemple, $\norm{\nabla f(\bar \w)}\leq \eps$ dans le cas sans contraintes). Notons que, en appentissage, avoir une précision $\eps$ très petite ne fait aucun sens, d'après les discussions de la section précédente.

Dernière généralité à ne pas oublier : il est très difficile, voire impossible, de résoudre -- même approximativement -- le problème général~\eqref{eq:opt} avec $f$ et $C$ quelconque. Il est essentiel de disposer de structures mathématiques pour envisager de résoudre les problèmes associés. Dans notre cas, deux propriétés importantes sont les propriétés de différentiabilité et de convexité, que nous allons revoir rapidement à la section suivante. 

%======================================================================
\paragraph{En théorie}

L'optimisation est une discipline des mathématiques appliquées qui étudie les algorithmes calculant des solutions approchées au problème \eqref{eq:opt}. On cherche à établir des garanties de convergence et des garanties d'efficacité numérique, ainsi que des garanties statistiques sur la solution obtenue. Cet aspect théorique n'est pas développé ici. Nous ne faisons que mentionner, en section \ref{sec:algo}, quelques résultats sur l'analyse de convergence des algorithmes, qui est une question centrale: il s'agit de montrer la convergence vers une solution optimale et aussi la vitesse de cette convergence, afin d'avoir des estimations de la complexité théorique des algorithmes dans certaines situations particulières. 


%======================================================================
\paragraph{En pratique}

Il existe de nombreux algorithmes et des logiciels spécialisés pour certaines classes de problèmes d'optimisation. Il existe aussi des logiciels \guil{modeleurs} qui font une interface entre les utilisateurs et les logiciels d'optimisation, permettant leur utilisation même par les non-experts. Par exemple, en science des données, la bibliothèque\;\texttt{scikit-learn} (en Python) rassemble des méthodes pour l'apprentissage utilisant un format unique qui permet de les utiliser avec une connaissance limitée des algorithmes d'optimisation sous-jacents. En général, il faut absolument se servir de ces algorithmes et de ces logiciels généralistes lors de la modélisation et pour résoudre les problèmes pour lesquels la performance numérique n'est pas la priorité. Pour les problèmes de grande taille ou dans un cadre exigeant, il est par contre fondamental de développer des algorithmes spécifiques exploitant la structure des problèmes visés. Dans tous les cas, il est utile de connaître les idées de base de l'optimisation, pour interpréter les sorties des logiciels ou régler les paramètres des algorithmes. C'est l'objectif de ce chapitre de présenter ces idées de base dans un cadre d'analyse\;de\;données.


%======================================================================
\subsection{Analyse convexe pour impatients}
%======================================================================

Les problèmes d'optimisation convexe forment une classe de problèmes que l'on peut bien étudier et résoudre, en théorie et en pratique. De plus, les méthodes pour le cas non convexe sont souvent des extensions de méthodes convexes ou utilisent des sous-problèmes convexes. Nous présentons ici, de manière aussi légère que possible, le minimum d'analyse convexe pour comprendre les algorithmes et les modèles en science des~données.

%======================================================================
\paragraph{Convexité et différentiabilité}

Une fonction $f\colon\RR^d\rightarrow \RR$ est dite convexe si elle vérifie pour tous $\w$ et $\v$:\index{convexité}
\[
f(\al \w + (1-\al)\v) \leq \al f(\w) + (1-\al)f(\v) \qquad \text{pour tout}~ \al\in [0,1]
\]
Cette inégalité signifie géométriquement que le graphe de $f$ est en-dessous du segment rejoignant 
$(\w,f(\w))$ et $(\v, f(\v))$ dans $\RR^d\times\RR$, comme l'illustre la figure\;\ref{fig:conv} avec $d=1$. 
Par exemple, $f(\om)= \exp(\om)$ et $f(\om)= \om\log(\om)$ (pour $\om>0$) sont des fonctions convexes d'une variable réelle,
de même que les fonctions pertes $\ell(y,\cdot)$ vues à la section\;\ref{sec:gene}. Les fonctions affines et les normes (comme $\norm{\cdot}_2$ et $\norm{\cdot}_1$) sont convexes en toute dimension. On montre souvent qu'une fonction $f$ est convexe sans revenir à la définition en justifiant que $f$ est construite à partir de fonctions convexes de base et d'opérations qui préservent la convexité (voir exercice 5.1). 

\begin{figure}[ht!]
  \begin{minipage}[c]{0.49\linewidth}
      \includegraphics[width=\linewidth]{Figs/chapter-Optimisation/convex.pdf}      
    \end{minipage} \hfill
    \begin{minipage}[c]{0.45\linewidth}
      \includegraphics[width=\linewidth]{Figs/chapter-Optimisation/not-convex.pdf}
    \end{minipage}
\caption{Illustration de la convexité avec des fonctions de $\RR$ dans $\RR$. \`A gauche, une fonction convexe; à droite, une fonction non convexe (dont on note qu'elle admet un minimum local en $\om$ qui n'est pas un minimum).\label{fig:conv}}
\end{figure}


La convexité entraîne des propriétés particulières sur les gradients (et les matrices hessiennes) des fonctions convexes différentiables
\footnote{Il existe toute une théorie qui étudie le cas des fonctions non différentiables, avec en particulier un objet qui explique le premier ordre des fonctions convexes (le sous-différentiel et les sous-gradients). Nous ne parlons pas de ces notions ici.}. % et référons au livre \cite{hiriart:2001} pour une présentation détaillée.}.
Par exemple, une propriété caractéristique des fonctions convexes est la suivante (illustrée sur la figure\;\ref{fig:tang}):
\begin{equation}\label{eq:tang}
f(\w) + \ip{\nabla f(\w)}{(\v-\w)} \leq f(\v) \qquad \text{pour tous~$\w$ et $\v$}
\end{equation}
Cette propriété implique, en particulier, qu'un minimum local ou même un point critique d'une fonction convexe est en fait un minimum global: on voit aisément que la propriété $\nabla f(\w)=0$ implique que $\w$ est un minimum global de~$f$. 
En termes mathématiques: la condition d'optimalité \eqref{eq:condopt}, toujours nécessaire, est aussi suffisante dans le cas convexe.

\begin{figure}[ht!]
\centering
\begin{minipage}[c]{0.48\linewidth}
      \includegraphics[width=\linewidth]{Figs/chapter-Optimisation/tangente.pdf}
\end{minipage}
\caption{Illustration de l'inégalité \eqref{eq:tang} pour une fonction convexe de $\RR$ dans~$\RR$. Une fonction convexe est toujours au-dessus de ses tangentes.\label{fig:tang}}
\end{figure}

%======================================================================
\paragraph{Convexité forte et différentiabilité forte}

Certaines fonctions convexes ont des propriétés supplémentaires utiles pour l'optimisation. Parmi ces dernières, la convexité forte mérite une attention particulière, car elle garantit l'existence et l'unicité d'une solution optimale, ainsi que de meilleures propriétés des algorithmes.\index{convexité!forte}
 
On dit qu'une fonction est fortement convexe si \guil{elle reste convexe quand on lui enlève un peu de convexité quadratique}; plus précisément, $f$ est fortement convexe (de module $\mu>0$) si la fonction $\w\mapsto f(\w) - \mu{\norm{\w}_2}^{\!2}$ est convexe. On démontre alors aisément qu'une fonction fortement convexe admet un unique minimum.
%Inversement on dit qu'un fonction est faiblement convexe si elle devient convexe quand on lui ajoute un peu de convexité":
%$f$ est faiblement convexe (de paramètre $\al>0$) si la fonction $x\mapsto f(x) + \al{\norm{x}_2}^{\!2}$ est convexe.
Il est fréquent en analyse de données que les problèmes d'optimisation soient construits pour être en fait fortement convexes, en choisissant lors de la modélisation d'ajouter d'une régularisation par la norme $\ell_2$. C'est par exemple le cas pour les moindres carrés avec régularisation $\ell_2$ de \eqref{eq:ridge} où la fonction est fortement convexe de module $\mu=\la$. Notons aussi qu'une fonction non différentiable peut aussi être fortement convexe, par exemple $R(\w) = \norm{\w}_1+\frac{\mu}{2}{\norm{\w}_2}^{\!2}$ appelée parfois la régularisation \guil{elastic-net}.

Une autre propriété très utile en optimisation est la différentiabilité forte de l'objectif, qui correspond intuitivement à ce que la fonction ne puisse pas \guil{changer de forme trop rapidement}. Plus précisément, on dit que $f$~est $L$-fortement différentiable (ou $L$-lisse) si elle est de classe $C^1$ avec un gradient Lipschitz de constante $L>0$, \cad :
\[
\norm{\nabla f(\w)- \nabla f(\v)}_2\leq L\norm{\w-\v}_2\qquad \text{pour tous $\w$ et $\v$ dans $\RR^d$}
\]
Par exemple, on montre aisément que la fonction-objectif de \eqref{eq:ridge} est fortement différentiable avec une constante $L$ inférieure à la somme de $\la$ et de la plus grande valeur propre de la matrice $\trans{X}X$. Notons que l'impact de ces deux notions\footnote{Mentionnons que les deux notions de convexité forte et différentiabilité forte sont reliées par une jolie \nom{dualité}. Nous disposons en effet du résultat suivant :
\[
\text{$f$ est $L$-fortement différentiable} \iff 
%si et seulement si 
\text{$f^*$ est $1/L$-fortement convexe}\]
mettant en jeu la fonction convexe conjuguée $f^*$ définie comme le résultat du problème d'optimisation suivant (paramétré par $\v$)
$
f^*(\v) = \max_{\w} ~\ip{\w}{\v} -f(\w).
$}
de convexité/différentiabilité forte porte sur les résultats mathématiques et aussi sur le comportement des algorithmes, comme nous allons le voir dans la section suivante. 

%Terminons en mentionnant que ces deux notions de convexité forte et différentiabilité forte sont reliés par une jolie \nom{dualité} en analyse convexe. Nous disposons en effet du résultat suivant 
%\[
%\text{$f$ est $L$-fortement différentiable} \iff 
%%si et seulement si 
%\text{$f^*$ est $1/L$-fortement convexe}\]
%mettant en jeu la fonction convexe conjuguée $f^*$ définie comme le résultat du problème d'optimisation suivant (paramétré par $\v$)
%\[
%f^*(\v) = \max_{\w} ~\ip{\w}{\v} -f(\w).
%\]
%% pour une fonction $f$ convexe continue,




%======================================================================
\subsection{Algorithmes d'optimisation}\label{sec:algo}
%======================================================================


Cette section donne une vue d'ensemble des méthodes du premier ordre en optimisation. Nous insistons sur la méthode de descente de gradient et ses généralisations pour les problèmes d'optimisation qui ont typiquement la forme suivante :
\begin{equation}\label{eq:pbcom}
\accol{\min \quad f(\w) + g(\w)\\ \quad \w\in \RR^d}
\end{equation}
où la fonction-objectif est une somme de deux fonctions de natures différentes, 
%\footnote{Ce problème inclut le cas avec contraintes \eqref{eq:opt} en prenant pour $g$ la fonction indicatrice des contraintes, qui vaut $g(\w)=0$ si $\w\in C$ et $g(\w)=+\infty$ si $\w\not\in C$.}
comme pour \eqref{eq:minERM}.
La présentation des algorithmes se déroule en trois parties traitant tour à tour le cas d'une seule fonction $f$ différentiable, puis d'une somme avec $f$ différentiable et $g$ qui ne l'est pas partout, et enfin d'une somme avec $f$ et $g$ fonctions non différentiables.


%======================================================================
\paragraph{Fonction-objectif différentiable: algorithme du gradient}

Commençons par le cas basique où $g=0$ et $f$ est différentiable. Géométriquement, le gradient $\nabla f(\w)$ en un point $\w$ donne la direction de la plus forte pente locale autour de $\w$. Il est donc naturel de suivre cette direction pour essayer de diminuer la valeur $f$; c'est ce que fait l'algorithme du gradient, dont l'itération s'écrit:\index{algorithme!du gradient}\index{descente|see{gradient}}
\begin{equation}\label{eq:grad}
\w_{k+1} = \w_k -\ga_k\nabla f(\w_k)
\end{equation}
Ici, $k$ est le compteur de l'itération et $\ga_k$ est un pas de descente, choisi pour garantir la convergence de l'algorithme.\index{pas de descente}\index{taux d'apprentissage} \index{pas d'apprentissage|see{pas de descente}}\index{learning rate|see{pas de descente}}
Le choix de ce pas de descente (appelé \nom{taux d'apprentissage} dans les applications en apprentissage) est fondamental pour l'efficacité de l'algorithme: pour des mauvais choix de pas, la méthode peut être extrêmement lente ou même ne pas converger (voir exercice 5.2). En pratique, le coût  de cette iteration en termes de temps de calcul se réduit essentiellement au coût de calcul du gradient.

Cet algorithme est suffisamment simple pour admettre une analyse mathématique complète. % dans le cas où $f$ est convexe.
%Avec quelques hypothèses sur $f$, 
En effet, il est possible d'analyser le nombre d'itérations nécessaires pour atteindre une solution approchée à $\eps$ près, pour toute fonction convexe $f$ satisfaisant certaines hypothèses. Une hypothèse standard, vérifiée 
dans presque toutes les applications en analyse de données, est que la fonction objectif $f$ soit fortement différentiable avec une constante $L$ (parfois calculée a priori et souvent estimée au fil des itérations). Sous cette hypothèse, on peut alors démontrer que la convergence est de l'ordre de $1/k$ en pire cas; plus précisément, pour toute fonction $f$ $L$-fortement différentiable et tout point initial $\w_0$, l'algorithme du gradient \eqref{eq:grad} avec un pas constant $\ga_k=1/L$ produit à l'itération $k$ un point $\w_k$ satisfaisant:\index{convergence}
\[
f(\w_k) -f_\star \leq \frac{2}{k}{\norm{\w_0-\w_\star}_2}^2
\]
Si la distance de l'itéré initial $\w_0$ à une solution optimale $\w_\star$ (qui apparaît ci-dessus) peut être bornée par une constante $D$, on obtient alors la garantie d'avoir une solution approchée à $\eps$ près, en $2D/\eps$ itérations (et ce, quelle que soit $f$).

Sans aucun réglage supplémentaire, cet algorithme s'avère même plus rapide si la fonction est de plus $\mu$-fortement convexe: l'algorithme du gradient (toujours avec le même pas $\ga_k=1/L$) converge \emph{de manière exponentielle} vers l'unique solution optimale:\index{convergence!exponentielle}
\begin{equation}\label{eq:cvg-expo}
{\norm{\w_k-\w_\star}_2}^2 \leq (1-\mu/L)^k {\norm{\w_0-\w_\star}_2}^2 \leq e^{-k\mu/L} {\norm{\w_0-\w_\star}_2}^2
\end{equation}
%Le terme de droite s'écrit comme $e^{-t/\kappa}$ (à une constante multiplicative près). 
Ce résultat théorique révèle le ratio $\kappa=\mu/L\leq 1$ qui joue le rôle de \nom{conditionnement}\index{conditionnement}
pour le problème de minimisation (d'une fonction $\mu$-convexe et $L$-lisse), comme on le voit sur le dessin de la figure~\ref{fig:ratio}. Dans le cas d'une fonction de régression linéaire $f(\w)={\norm{X\w-\y}_2}^2$, on montre que $\kappa$ correspond au ratio des valeurs propres extrêmes $\kappa=\lambda_{\min}(\trans{X}X)/\lambda_{\max}(\trans{X}X)$.
\begin{figure}[h!]
\centering
    \begin{minipage}[c]{0.85\linewidth}
      \includegraphics[width=\linewidth]{Figs/chapter-Optimisation/conditionnement.pdf}
    \end{minipage}
    \caption{Dessin des lignes de niveau de deux fonctions de $\RR^2$ dans $\RR$ avec les itérés de l'algorithme du gradient, illustrant la vitesse de convergence \eqref{eq:cvg-expo} et le rôle du conditionnement $\kappa=\mu/L$.
    Quand $\kappa$ est grand (comme à droite), le problème est bien conditionné et l'algorithme converge vite; plus le ratio est petit (comme à gauche), plus le problème est dur et plus l'algorithme patine.\label{fig:ratio}}
\end{figure}

%======================================================================
\paragraph{Ingrédients de base et recettes avancées}

L'algorithme de gradient est simple à comprendre et à coder. C'est en fait l'archétype des algorithmes d'optimisation, avec deux ingrédients essentiels: le pas de descente $\ga_k$ et la direction de descente donnée ici par le gradient. Chacun de ces deux ingrédients peut être amélioré:
\begin{itemize}
\item L'algorithme converge plus vite en choisissant des pas de descente particuliers pour une fonction ou une classe de fonctions. Par exemple, on peut améliorer la vitesse \eqref{eq:cvg-expo} dans le cas fortement convexe en choisissant un pas dépendant de la constante de convexité forte\footnote{En prenant un pas constant $\ga_k=2/(L+\mu)$ dépendant à la fois de $L$ la constante de différentiabilité forte et de $\mu$ la constante de convexité forte, la vitesse exponentielle de l'algorithme du gradient peut être améliorée en:
\begin{equation}\label{eq:cvg-expo2}
{\norm{\w_k-\w_\star}_2}^2 \leq \big(\frac{\kappa-1}{\kappa+1}\big)^k {\norm{\w_0-\w_\star}_2}^2 \leq \exp(-4k/(\kappa+1))) {\norm{\w_0-\w_\star}_2}^2
\end{equation}
Le gain est simplement de remplacer $\kappa$ par $(\kappa+1)/4$, ce qui n'est pas spectaculaire, mais permet de retrouver dans le cas particulier de $f$ quadratique un résultat bien connu en analyse numérique matricielle.
}. En général, c'est encore mieux en calculant des $\ga_k$ adaptatifs, spécifiques à une fonction et une itération particulière.  

\item D'autres directions de descente sont meilleures et permettent d'éviter l'effet de zig-zag que l'on voit sur la figure\;\ref{fig:ratio}. L'idée est d'incorporer dans l'algorithme de l'information supplémentaire sur le second ordre de $f$, comme dans la méthode de gradient accéléré par inertie\footnote{Les méthodes d'accélération par inertie ou relaxation ont fait l'objet de nombreuses études récentes dans le cadre de l'analyse de données. Une jolie introduction est proposée sur le site web \texttt{https://distill.pub/2017/momentum/}.} ou la méthode de Newton.\index{méthode!de Newton} % FORMULES ?
\end{itemize}

En particulier, les méthodes de quasi-Newton (BFGS) à mémoire limitée, combinant les deux aspects (pas adaptatif et information supplémentaire), sont reconnues comme efficaces et sont utilisées depuis longtemps pour des problèmes avec données météorologiques.  En général, les algorithmes du second ordre (comme des méthodes de type Newton), incorporant de l'information sur la courbure de la fonction, demandent moins d'itérations pour atteindre une bonne précision. En contrepartie, ils nécessitent souvent de connaître et de manipuler cette information supplémentaire (avec notamment des calculs plus coûteux à chaque itération), ce qui n'est rentable que si l'on exploite la structure particulière des problèmes d'optimisation en science des données (de grande taille et n'exigeant pas une précision fine).
 

%
%\[
%\w_{k+1} = \w_k -\ga_kW_k\nabla f(\w_k)
%\]
%Avec $W_k$ une matrice qui approche l'inverse de $\nabla^2 f(\w_k)$ dans les bonnes directions.
%
%pour obtenir une precision à $\eps$ pres
%
%gradient $O(\frac{L}{\mu}\log(1/\eps))$
%
%Nesterov $O(\sqrt{\frac{L}{\mu}}\log(1/\eps))$

%======================================================================
\paragraph{Fonction-objectif composée: algorithme de gradient-proximal}

Considérons maintenant le problème \eqref{eq:pbcom} où $f$ est différentiable et $g$ est une fonction convexe non différentiable simple. Dans ce cas, le gradient de la fonction-objectif n'existe plus partout et une application directe de l'algorithme du gradient \eqref{eq:grad} n'est pas possible (penser par exemple à la minimisation de la valeur absolue dans $\RR$). Le fait que $g$ soit \guil{simple} permet de quand même définir un algorithme de premier ordre qui exploite le fait que l'objectif est une somme.

%Les fonctions $g$ qui apparaissent dans les applications en science des données admettent souvent la propriété particulière suivante.  
On dit que $g$ est \nom{prox-simple} si on peut aisément calculer (via une expression explicite peu co\^uteuse ou une sous-routine très efficace) la solution optimale du problème d'optimisation suivant, pour un point $\w$ fixé:\index{opérateur proximal}\index{prox-simple|see{opérateur proximal}}
\begin{equation}\label{eq:prox}
\min_{\z\in \RR^d} ~~g(\z) +\frac{1}{2}\norm{\z-\w}^2
\end{equation}
Les exemples en apprentissage de telle fonctions sont légions: régularisation par norme $\ell_1$ (voir exercice 5.3), norme $\ell_1$ par blocs, contraintes de type bo{\^{\i}}te... 
%Nous dévelopons l'exemple de la norme $\ell_1$ ci-dessous. 
Observons sur \eqref{eq:prox} que le terme quadratique additionnel garantit qu'il existe une unique solution pour tout paramètre $\w$; cela permet de définir l'opérateur proximal\index{opérateur proximal} de $g$ noté $\prox_g\colon \RR^d\rightarrow\RR^d$ qui à $\w$ associe cette unique solution de \eqref{eq:prox}.

%\begin{exem}{La régularisation $\norm{\cdot}_1$ est prox-simple}
%Traitons explicitement l'exemple de 
%D'autres régularisations construite à partie de la norme $\ell_1$ admettent aussi des ope
%Norme matricielle comme la "norme-trace" ou norme nucléaire $\norm{X}_*=\norm{\si(X)}_1$ la norme $\ell_1$ des valeurs singulières de la matrice.
%\end{exem}
%

L'algorithme du gradient-proximal consiste à faire un pas de gradient sur~$f$ suivi d'une correction par l'opérateur proximal\index{opérateur proximal} associé à $g$, ce qui donne l'itération suivante:\index{algorithme!de gradient proximal}
\begin{equation}\label{eq:proxgrad}
\w_{k+1} = \prox_{\ga_k\,g}(\w_k -\ga_k \nabla f(\w_k))
\end{equation}
Notons que, pour $g=0$, l'opérateur proximal\index{opérateur proximal} est l'identité et on retombe sur l'algorithme du gradient discuté précédemment. 
Il se trouve que l'analyse de convergence de cet algorithme (ainsi que de ses extensions et accélérations) donne des résultats très similaires à ceux de l'algorithme du gradient rappelés juste avant.

Plaçons-nous finalement dans le cas général de \eqref{eq:pbcom} quand $f$ et $g$ sont toutes deux non différentiables et prox-simples. Comme l'opérateur proximal\index{opérateur proximal} d'une somme n'est pas simple en général, nous devons reformuler le problème pour exploiter chaque opérateur proximal indépendamment. La manipulation consiste à ajouter une variable supplémentaire $z\in \RR^d$ pour reformuler \eqref{eq:pbcom} comme :
\[
\accol{\min \quad f(\w) + g(\v)\\ \quad \w=\v}%_{\w\in \RR^d,z\in \RR^d}
\]
Intuitivement, on peut alors faire une itération proximale sur chaque variable en ajoutant une étape qui connecte les deux, ce qui donne:\index{ADMM}
\begin{equation}\label{eq:admm}
\accol{\w_{k+1} = \prox_{f/\rho} (\v_k - \u_k/\rho)\\
\v_{k+1} = \prox_{g/\rho} (\w_k + \u_k/\rho)\\
\u_{k+1}=\u_k + \rho(\w_{k+1}-\v_{k+1})}
\end{equation}
pour $\rho>0$ quelconque fixé.
Cette itération correspond à celle d'un algorithme appelé ADMM (pour \nom{Alternating Direction Method of Multipliers}\footnote{La variable $u$ s'interprète en effet comme une variable duale. Une justification propre de cet algorithme utilise la dualité de l'optimisation convexe.}) dans un cas particulier. 
%Notons simplement ici qu'aucune hypothèse particulière n'est nécessaire pour la convergence de cet algorithme.

Terminons en disant que, si les fonctions non différentiables apparaissant dans \eqref{eq:pbcom} ne sont pas prox-simples, d'autres familles d'algorithmes existent, parmi lesquelles les plus populaires en analyse de données sont :
\begin{itemize}
\item les méthodes de gradient conditionnel (appelé aussi Frank-Wolfe);
\item les méthodes proximales dans des géométries favorables (utilisant la notion de \nom{divergence de Bregman} pour mesurer l'éloignement).
\end{itemize}


%======================================================================
\paragraph{Bilan sur les algorithmes d'optimisation de cette section}

Les algorithmes présentés dans cette section sont simples et permettent de présenter de nombreux aspects de l'optimisation (descente, pas, propriétés de convergence, rôle de la convexité, de la différentiabilité). Ils ont aussi de bonnes propriétés pour les passages à l'échelle présentées dans la section suivante. Notons enfin qu'ils se trouvent remarquablement robustes à l'utilisation d'approximation stochastique des gradients et des opérateurs proximaux,\index{opérateur proximal} et s'accommodent des techniques de randomisation.% (par moyennage, échantillonnage, ou tirage de Monte-Carlo...).

%XX  utilisées par exemples pour les systèmes de recommendation utilisant la norme-trace




%======================================================================
%======================================================================
\section{Algorithmes en science des données}\label{sec:appli}
%======================================================================
%======================================================================

Comme vu à la section \ref{sec:gene}, de nombreux problèmes d'optimisation en science des données s'écrivent
comme la minimisation de sommes de fonctions du type :% \eqref{eq:eq:sumfinie}
\begin{equation}\label{eq:sumfinie2}
\min_{\w\in \RR^d} \quad \frac{1}{n}\sum^n_{i=1} f_i(\w) + r(\w)
\end{equation}
Les fonctions $f_i\colon\RR^d\rightarrow\RR$ sont similaires et concernent chacune une partie des données: typiquement, $f_i$ mesure la qualité de la prédiction associée à $\w$ sur un bloc de données. 
Pour fixer les idées, on peut garder à l'esprit  la version basique $f_i(\w) = \ell\big((\x_i,y_i),\w\big)$,
pour laquelle \eqref{eq:sumfinie} correspond à \eqref{eq:minERM}. 
Dans le cas de blocs de données (par exemple lorsque les données sont stockées sur différentes machines), on se trouve alors avec des fonctions du type:
\[
f_i(\w) = \sum_{j\in B_i}\ell\big(y_i,h(\x_j,\w)\big),\quad \text{où $B_i$ désigne un bloc de données.}
\]

Ce problème d'optimisation \eqref{eq:sumfinie2} pourrait être vu comme étant de la forme \eqref{eq:pbcom} avec:
\begin{equation}\label{eq:fg}
\text{$f(\w) = \frac{1}{n}\sum^n_{i=1} f_i(\w)$ ~~~et~~~ $g(\w) = r(\w)$}
\end{equation}
et ainsi pourrait être résolu par les algorithmes de la section\;\ref{sec:algo}.
Malheureusement, une utilisation directe de ces derniers n'est pas appropriée dans le cas de données volumineuses, pour lesquelles certaines opérations élémentaires, demandant l'utilisation de toutes les données, s'avèrent très coûteuses.

Illustrons ceci dans le cas où les $f_i$ sont différentiables et $r$ prox-simple, comme pour \eqref{eq:lasso}. 
%Regardons par exemple le problème de régression linéaire parcimonieuse \eqref{eq:lasso}. 
%Dans ce cas, $f$ est différentiable et $g$ est prox-simple: 
On peut utiliser l'algorithme du gradient proximal \eqref{eq:proxgrad} pour résoudre le problème \eqref{eq:sumfinie2}, ce qui donne l'itération:
\begin{equation}\label{eq:gradproxfini}
\w_{k+1} = \prox_{\ga_k r}\left(\w_k -\frac{\ga_k}{n}\sum^n_{i=1}\nabla f_i(\w_k)\right)
\end{equation}
Pour de nombreuses régularisations classiques, le calcul de l'opérateur proximal\index{opérateur proximal} n'est pas trop coûteux, ni en temps de calcul, ni en mémoire (voir par exemple $\prox_{\norm{\cdot}_1}$ en \eqref{eq:soft}). La principale opération de cette itération est le calcul du gradient de la somme, qui requiert une passe sur toutes les $n$ données pour calculer les gradients $\nabla f_i(\w_k)$, auquel s'ajoute le coût de sommer $n$ vecteurs de $\RR^d$. C'est cette opération qui est coûteuse quand $n$ et $d$ sont grands et qui rend cet algorithme peu efficace.

L'idée pour passer à l'échelle sur des données volumineuses est de traiter chaque bloc de données indépendamment, plutôt que de considérer toutes les données en même temps. Nous présentons deux stratégies allant dans ce sens: les algorithmes incrémentaux, pour lesquels une itération correspond à l'exploitation d'un seul bloc de données (section \ref{sec:sto}),
et les algorithmes distribués, adaptés au cas où les données sont stockées sur des machines différentes (section \ref{sec:distrib}).
%\footnote{Nous séparons les deux types pour en présenter les idées sous-jacentes; dans les applications pratiques, il y a souvent un mélange des deux.} 



%======================================================================
\subsection{Algorithmes incrémentaux}\label{sec:sto}
%======================================================================

Les algorithmes incrémentaux opèrent une mise à jour à chaque lecture d'un bloc de données (au lieu de tous les utiliser en même temps). Nous présentons le plus simple d'entre eux, l'algorithme du gradient incrémental (ou gradient stochastique), puis nous discutons ses extensions.

%======================================================================
\paragraph{Une donnée, une itération}

Plaçons-nous dans cette section dans le cas basique de la minimisation de fonctions différentiables ($r=0$ dans \eqref{eq:sumfinie2}). %Par exemple, logistique régularisée par norme 2. XXX
Dans ce cas, on pourrait utiliser l'algorithme du gradient :
\begin{equation}\label{eq:gradfini}
\w_{k+1} = \w_k -\frac{\ga_k}{n}\sum^n_{i=1}\nabla f_i(\w_k)
\end{equation}
Comme expliqué précédemment pour \eqref{eq:gradproxfini}, l'inconvénient est qu'il nécessite de 
faire, à chaque itération, un passage sur toutes les données afin de calculer tous les $\nabla f_i(\w_k)$.
L'idée est alors de faire une itération en utilisant un seul bloc de données et ainsi un seul $\nabla f_i(\w_k)$. 

\`A l'itération $k$, on choisit un indice $i_k$, puis on calcule $\nabla f_{i_k}(\w_k)$ sur le bloc de données correspondant à $i_k$ et on procède à la mise à jour:\index{algorithme!du gradient stochastique} 
\begin{equation}\label{eq:gradsto}
\w_{k+1} = \w_k -\ga_k\nabla f_{i_k}(\w_k) \qquad\text{$i_k\in \{1,\ldots,n\}$}
\end{equation}
Cette itération \eqref{eq:gradsto} correspond à l'algorithme du gradient incrémental, dépendant 
du choix de $i_k$ (qui peut se faire de manière cyclique, gloutonne ou aléatoire).\index{algorithme!aléatoire}
%Si on suppose que les données sont i.i.d, on obtient une meilleure convergence si $i_k$ est choisi de manière aléatoire
En exploitant la redondance\footnote{Considérons une situation dégénérée pour bien comprendre l'idée. Supposons que l'ensemble de données est constitué de dix copies identiques du même bloc; autrement dit, qu'il y a dix $f_i$ tous égaux. Dans ce cas, l'algorithme~\eqref{eq:gradsto} donne, pour un coût de calcul $n$ fois plus petit, exactement la même itération que l'algorithme du gradient~\eqref{eq:gradfini}.  \'Evidemment, en pratique, les données ne sont pas dupliquées comme dans cette situation. Néanmoins, beaucoup de problèmes d'analyse de données contiennent une bonne partie de données redondantes, ce qui rend cette approche intéressante.} inhérente des données, cet algorithme donne une itération de minimisation pour un coût de calcul nettement moindre que \eqref{eq:gradfini}: le calcul de $\nabla f_{i_k}(\w_k)$ est bien moins coûteux que celui de $\nabla f(\w)$ tout entier (typiquement, $n$~fois moins coûteux).\index{données redondantes}


%======================================================================
\paragraph{Optimisation par gradient stochastique}


Dans le cas où $i_k$ est aléatoire et choisi uniformément dans $\{1,\ldots,n\}$, l'itération \eqref{eq:gradfini} s'appelle aussi gradient stochastique: on observe en effet que l'espérance de la direction de descente (aléatoire) correspond au gradient :
\begin{equation}\label{eq:esp}
\mathbb{E}[\nabla f_{i_k}(\w_k)] =\nabla f(\w_k)
\end{equation}
%Plus généralement
%\[
%\w_{k+1} = \w_k -\ga_k g_k 
%\]
%avec $g_k\in \RR^n$ une direction de descente aléatoire telle que $\mathbb{E}[g_k] =\nabla f(\w_k)$
Utiliser une direction de descente aléatoire amène une décroissance aléatoire (voir la figure\;\ref{fig:gradsto}). On peut tout de même contrôler la convergence un prenant un pas $\ga_k$ qui tend vers $0$ (mais pas trop vite: typiquement $\ga_k=1/k$).
Prendre des pas de plus en plus petits entra{\^{\i}}ne néanmoins que cet algorithme converge lentement.  %On démontre alors la convergence (de $\mathbb{E}[\norm{\nabla f(\w_k)}_2}^2]$) vers un point $\bar \w$ tel que $\nabla f(\bar \w) =0 $ (qui est un minimum si $f$ est convexe). 
%XXX DESSIN.
\begin{figure}[ht!]
\centering
 \begin{minipage}[c]{0.49\linewidth}
      \includegraphics[width=\linewidth]{Figs/chapter-Optimisation/gradient.pdf}      
    \end{minipage} \hfill
    \begin{minipage}[c]{0.48\linewidth}
      \includegraphics[width=\linewidth]{Figs/chapter-Optimisation/SG.pdf}
    \end{minipage}
    \caption{Itérations des algorithmes du gradient (à gauche) et du gradient stochastique (à droite) sur les lignes de niveau d'une fonction de $\RR^2$ dans $\RR$. L'espérance des directions de descente correspond à la direction du gradient et on récupère des garanties de convergence si on peut aussi contrôler la variance des directions. \label{fig:gradsto}}
\end{figure}

Par rapport à l'algorithme du gradient, on a ainsi dans une situation opposée 
concernant l'équilibre \guil{coût de calcul} vs \guil{vitesse de convergence}:

\begin{itemize}
\item pour l'algorithme du gradient \eqref{eq:gradfini}:  un coût d'itération proportionel à $n$, mais une vitesse de convergence de l'ordre de $e^{−k/\kappa}$ dans le cas fortement convexe (voir \eqref{eq:cvg-expo});
\item pour l'algorithme du gradient stochastique avec pas $\ga_k=1/k$: un coût d'itération indépendant de $n$, mais une vitesse de convergence de l'ordre  $\kappa/k$ uniquement.
\end{itemize}

Ainsi, même si chaque itération est bien moins coûteuse que pour la méthode du gradient, l'algorithme du gradient stochastique n'est pas efficace en fin de convergence. Son intérêt réside dans ses premières itérations: on observe en effet systématiquement une décroissance initiale rapide. En pratique, il est encore une fois fondamental de régler correctement\footnote{Deux recettes à garder en tête : prendre un pas initial modéré par la constante de Lipschitz pour garantir une décroissance dès les premières itérations; ensuite, prendre des pas aussi grands que possible (dans les limites de résultats généraux de convergence ou même au-delà si on a besoin de faire décro{\^{\i}}tre rapidement l'objectif).} le pas de descente $\ga_k$ pour observer une bonne convergence empirique. Au-delà de ces réglages empiriques, d'autres techniques de réduction de variance ont été proposées pour améliorer les performances de l'algorithme de gradient stochastique.

%mini-batch
%\[
%\w_{k+1} = \w_k -\frac{\ga_k}{|B_k|}\sum_{i\in B_k}\nabla f_i(\w)
%\]
%avec un sous-ensemble des données 
%blocs $B_k\in \{1,\ldots,n\}$
%et $|B_k|$ est la taille du bloc $k$.
%


%======================================================================
\paragraph{Réduction de l'aléatoire par moyennes}

Il existe des alternatives à l'utilisation de pas décroissants pour contrôler la variance\index{réduction de variance} dans \eqref{eq:esp} par moyennage\footnote{Il existe aussi des méthodes de réduction de variance utilisant la dualité.}.
Une première technique classique est de considérer les moyennes des itérés successifs produits par \eqref{eq:gradsto}.
Ce moyennage :
\[
\bar \w_k = \frac{1}{k} \sum^k_{\ell=1} \w_\ell
\]
amène une certaine stabilité de la suite $(\bar \w_k)_k$ et permet ainsi de prendre des pas plus grands.
Une seconde technique plus récente est de faire des moyennes des gradients successifs.
Il s'agit de garder en mémoire certains gradients calculés aux itérations précédentes et de les combiner pour obtenir une meilleure descente. Par exemple, on peut garder en mémoire un gradient par fonction $f_i$: à l'itération $k$, on tire aléatoirement l'indice $i_k$, on calcule 
$\nabla f_{i_k}(\w_k)$ et on met à jour la mémoire des $n$ gradients stockés:
\[
g_i^k = \nabla f_i(\w_k) ~~\text{si $i=i_k$ ~~et~~ $g_i^k = g_i^{k-1}$ sinon.}
\]
L'itération de l'algorithme appelé SAG (\textit{stochastic averaged gradient}) consiste alors à prendre la moyenne des gradients en mémoire:\index{algorithme!SAG}
\[
\w_{k+1} = \w_k -\frac{\ga_k}{n}\sum^n_{i=1}g_i^k
\]
On voit que cette itération peut se réécrire de la manière suivante, en remplaçant le vecteur gradient $i_k$ dans la moyenne des $g^k_i$\,:
\[
\w_{k+1} = \w_k -\ga_k\Big(\frac{1}{n}\sum^n_{i=1}g_i^k + \frac{1}{n}(\nabla f_{i_k}(\w_k) - g_{i_k}^k)\Big)
\]
De simples modifications de cette itération amènent à des algorithmes efficaces, ayant de plus de bonnes garanties théoriques. Par exemple\footnote{Il existe d'autres méthodes par moyennage, dont une, appelée SVRG, qui admet une implémentation avec faible mémoire, allégeant ainsi le coût supplémentaire en mémoire de $n$ vecteurs de $\RR^d$ des méthodes présentées ici.}, en supprimant le $1/n$ (ce qui s'interprète comme insister sur l'itération courante), on obtient l'itération de l'algorithme appelé SAGA\index{algorithme!SAGA}
\begin{eqnarray}
\w_{k+1} &=& \w_k -\ga_k\Big(\frac{1}{n}\sum^n_{i=1}g_i^k + (\nabla f_{i_k}(\w_k) - g_{i_k}^k)\Big) \nonumber \\
&=& \w_k -\ga_k\nabla f_{i_k}(\w_k) + \Big(\frac{1}{n}\sum^n_{i=1}g_i^k  - g_{i_k}^k\Big)\label{eq:saga}
\end{eqnarray}
%\[
%\w_{k+1} = \w_k -\ga_k\nabla f_{i_k}(\w_k) + \Big(\frac{1}{n}\sum^n_{i=1}g_i^k  - g_{i_k}^k\Big)
%\]
Avec une meilleure direction de descente stochastique que celle de \eqref{eq:gradsto}, 
cet algorithme admet des pas de descente plus grands et s'avère ainsi plus rapide. %que le gradient stochastique en pratique
Dans le cas $\mu$-fortement convexe et $L$-fortement différentiable, on montre que \eqref{eq:saga} avec un pas constant $\ga_k=\frac{1}{2(\mu d + L)}$ admet la convergence exponentielle suivante (en erreur quadratique moyenne), qui très proche de ~\eqref{eq:cvg-expo2} pour l'algorithme du gradient:
\[
\mathbb{E}\left[{\norm{\w_k-\w_\star}_2}^2\right]\leq C \exp\left(\frac{-k\mu}{2(\mu d + L)}\right) %\qquad \text{(avec $C$ une constante dépendant de $\w_0$ et $\w_\star$)}
\]
avec une constante $C$ dépendant de $\w_0$ et $\w_\star$.

%méthode stochastique SG SAG SAGA SVRG 
%plus rapide que les méthode batch mais demande quand meme une mise a jour de $\w$ tout entier.



%======================================================================
\paragraph{Bilan sur les algorithmes incrémentaux}

Les méthodes du type gradient stochastique sont efficaces lors des premières itérations pour faire décro{\^{\i}}tre l'objectif rapidement sans avoir à exploiter la totalité des données, mais elles ralentissent très vite à cause d'un taux d'apprentissage décroissant imposé par la théorie. En acceptant un coût supplémentaire en mémoire et en calcul, elles sont modifiables pour accepter des pas constants et ainsi afficher une meilleure convergence en pratique et en théorie. 

Ces méthodes marchent particulièrement bien pour minimiser des fonctions fortement convexes, ce qui incite à ajouter, si nécessaire, une régularisation quadratique dans les modèles d'apprentissage, ce qui revient à prendre $f_i(\w) = \ell\big((\x_i,y_i),\w\big) + \la_2{\norm{\w}_2}^2$. Dans tous les cas, le réglage empirique des différents paramètres des algorithmes (notamment le pas de descente) requiert une attention particulière car il affecte fortement le comportement observé des algorithmes. 

\`A retenir : l'utilisation de l'aléatoire\footnote{D'autres méthodes aléatoires sont utiles en optimisation en grande dimension: par exemple, la mise à jour d'un bloc  aléatoire de variables, l'approximation des gradients ou des opérateurs proximaux\index{opérateur proximal} par des estimations stochastiques, ou l'accélération des routines de calcul matriciel par randomisation.} dans ces algorithmes incrémentaux permet de récupérer de l'information pertinente en ne traitant qu'une partie des données. Néanmoins, cet aléa requiert de pouvoir accéder rapidement à n'importe quelle partie des données. Ceci implique des restrictions sur les systèmes informatiques: les différentes machines de calcul doivent avoir une mémoire partagée contenant les données (en plus de ressources de calcul partagées). Si ce n'est pas le cas, d'autres algorithmes d'optimisation (déterministes) existent; ils sont présentés à la section suivante.

%Notons finalement qu'une itération des algorithmes incrémentaux, comme l'algorithme du gradient stochastique \eqref{eq:gradsto}, 

%Les méthodes de type gradient sont simples et permettent te traiter les problèmes de grande taille. Les méthodes de type gradient stochastique permettent de passer à l'échelle pour un grand nombre de données, avec malheureusement une convergence plus lente. Des extensions (comme SAG ou SVRG) permettent d'en gommer les défauts comme par exemple SAG pour le toute convergence, et etSVRG pour le problème de mémoire. Ses idées doit être vu comme des briques de base qui peuvent être étendues de nombreux problèmes plus complexes (non-lisse ou avec contraintes)


%======================================================================
\subsection{Algorithmes distribués}\label{sec:distrib}
%======================================================================

Dans cette section, nous considérons la situation où les données ne peuvent être stockées au même endroit, par exemple pour des raisons techniques (données trop volumineuses) ou des raisons juridiques (données privées, par exemple celles des téléphones portables). Dans ce cas, les données sont stockées sur différentes machines et la résolution d'un problème d'apprentissage \eqref{eq:sumfinie2} doit se faire avec des algorithmes d'optimisation \emph{distribuée}. Les questions reliées aux communications entre machines deviennent alors prépondérantes pour l'efficacité de ces algorithmes.\index{optimisation!distribuée}

Pour simplifier, plaçons-nous dans le cas où l'on dispose d'une architecture de calcul centralisée avec une machine-maître qui dirige $m$~machines-esclaves; pour fixer les idées, considérons que $m=n$,
%qu'il y a $n$ machines, \cad a autant de machines de blocs de données et de fonctions $f_i$. Les idées présente s'étendent 
\cad qu'il y a autant de machines que de blocs de données (et de fonctions $f_i$). 
L'indice $i$ réfère maintenant à la fois au bloc de données $i$ comme précédemment, mais aussi à la machine qui le stocke.

Il se trouve que les algorithmes du premier ordre rappelés en section~\ref{sec:algo} s'étendent naturellement dans ce cadre distribué, où les calculs coûteux sont gérés en local sur chaque machine, puis les résultats sont reliés entre eux par une étape de consensus. Cette section présente brièvement les idées de ces algorithmes distribués, d'abord dans le cas où les $f_i$ sont différentiables, puis dans le cas o\`u elles ne le sont pas partout.



%
%Comme présenté dans plusieurs chapitre de ce livre, les nouvelles technologies big data (logicielles et matérielles) connaissent de grandes avancées.
%
%exploiter les architectures de calculs complexes et les outils logiciels big data


%======================================================================
\paragraph{Calculer le gradient en parallèle}

Reprenons le problème d'optimisation \eqref{eq:sumfinie2}, qui est de la forme \eqref{eq:pbcom} avec les fonctions\;\eqref{eq:fg}. Supposons que les fonctions d'attache aux données $f_i$ sont différentiables, ce qui implique que $f$ l'est aussi avec :
\begin{equation}\label{eq:gradpara}
\nabla f(\w) = \frac{1}{n}\sum^n_{i=1}\nabla f_i(\w)
\end{equation}
On peut alors utiliser l'algorithme du gradient proximal qui donne l'itération \eqref{eq:gradproxfini}, dont l'opération coûteuse est le calcul de la somme des gradients.
Dans le cadre distribué de cette section, cette itération se prête bien au calcul parallèle: 
\begin{enumerate}
\item la machine-maître diffuse la variable $\w_k$ aux machines-esclaves; 
\item le calcul des $n$ gradients $\nabla f_i(\w_k)$ dépendant des données se fait en parallèle sur ces $n$ machines; 
\item la machine-maître centralise les résultats des calculs pour calculer $\w_{k+1}$, en faisant la moyenne des gradients, suivie de l'opérateur proximal\index{opérateur proximal} de~$r$.
\end{enumerate}
Ce schéma de calcul entre dans le cadre d'un formalisme \nom{map-reduce} (affectation-réduction) où la phase de \guil{map} distribue le calcul coûteux sur toutes les machines et la phase de réduction synchronise les résultats des calculs avec uniquement des opérations simples (figure~\ref{fig:mapreduce}).\index{map-reduce}
 
\begin{figure}[ht!]
\centering
\begin{minipage}[c]{0.75\linewidth}
      \includegraphics[width=\linewidth]{Figs/chapter-Optimisation/mapreduce.pdf}
\end{minipage}
\caption{Représentation schématique du cadre d'optimisation distribuée: une machine-maître coordonne les calculs effectués en parallèle sur les machines-esclaves, stockant chacune une partie des données. Soulignons que les algorithmes ne communiquent que des résultats de calcul, pas des données brutes. 
\label{fig:mapreduce}}
\end{figure}

Ce raisonnement, présenté ici dans un cadre centralisé, s'étend  à des systèmes informatiques plus complexes.
Dans le cas général, des logiciels \guil{big~data} permettent d'implementer cet algorithme sans trop se soucier du système sous-jacent. Par exemple, \texttt{Spark} permet de programmer de la même manière le cas $m<n$ que le cas $n=m$ considéré~ici.
 

 
%======================================================================
\paragraph{Calculer l'opérateur proximal en parallèle (après reformulation)}

\index{opérateur proximal}Toujours dans le cadre distribué avec le problème d'optimisation \eqref{eq:sumfinie2} (de la forme \eqref{eq:pbcom} avec \eqref{eq:fg}), supposons à présent que les fonctions\;$f_i$ sont non différentiables et prox-simples.
 % comme pour f_i(\w) = \norm{A_i\w-b_i}_1
Contrairement à la situation précédente, on n'est pas directement dans une situation de la section~\ref{sec:algo}: en effet, 
$f$ n'est pas prox-simple, car l'opérateur proximal d'une somme ne correspond pas à la somme des opérateurs proximaux\index{opérateur proximal} (contrairement aux gradients\;\eqref{eq:gradpara} dans le différentiable). On ne peut donc pas distribuer les calculs en appliquant directement l'algorithme ADMM \eqref{eq:admm} adapté à la minimisation de la somme de deux fonctions prox-simples. Une manière de contourner cette interdiction consiste à reformuler \eqref{eq:sumfinie2}, comme suit, sous une forme distribuable impliquant des fonctions prox-simples.
 
En donnant à chaque $f_i$ une copie locale de $\w$, notée $\w_i$,  le problème \eqref{eq:sumfinie2} peut s'écrire de manière équivalente comme :
 \begin{equation}\label{eq:formdistri}
\accol{\displaystyle\min_{(\w_1,\ldots,\w_n)\in \RR^{d\times n}} %~\frac{1}{n}\sum^n_{i=1} f_i(\w_i) +  r(\w)
~~=~~ \frac{1}{n}\sum^n_{i=1} \big(f_i(\w_i) +  r(\w_i) \big)\\
\qquad \w_1=\cdots= \w_n}
\end{equation}
En introduisant la fonction indicatrice de consensus $\delta\colon\RR^{d\times n}\rightarrow\RR\cup\{+\infty\}$ (définie par $\delta\big((\w_1,\ldots,\w_n)\big) = 0$ si $\w_1=\cdots= \w_n$ et $+\infty$ sinon), on voit que ce problème est de la forme \eqref{eq:pbcom}  avec :
\begin{align}
f\big((\w_1,\ldots,\w_n)\big) &= \frac{1}{n}\sum^n_{i=1}  f_i(\w_i) \label{eq:feclat}\\
%\quad\text{et}\quad 
g\big((\w_1,\ldots,\w_n)\big) &= \frac{1}{n}\sum^n_{i=1}  r(\w_i) + \delta\big((\w_1,\ldots,\w_n)\big) \label{eq:geclat}%\delta_{\w_1=\cdots= \w_n}\big((\w_1,\ldots,\w_n)\big).
\end{align}
On montre (voir exercice 5.4) que ces deux fonctions sont 
prox-simples et que leurs opérateurs proximaux\index{opérateur proximal} se décomposent selon chaque $\w_i$ et s'écrivent uniquement à l'aide des $\prox_{f_i}$ et de $\prox_{r}$.
On peut alors appliquer l'algorithme ADMM \eqref{eq:admm} avec les bonnes fonctions, ce qui donne après ré-écriture:\index{ADMM}
%des itérations éclatées sur chaque machine:
%pour toute machine $i$, 
\[
\accol{(\w_i)_{k+1} = \prox_{f_i/n\rho} (\v_k - (\u_i)_k/\rho) \qquad \text{pour toute machine $i$}\\
\v_{k+1} = \prox_{r/n\rho} \left(\frac{1}{n}\sum^n_{i=1}((\w_i)_{k+1} + (\u_i)_k/\rho)\right)\\
(\u_i)_{k+1}=(\u_i)_k + \rho((\w_i)_{k+1}-\v_{k+1}) \qquad \text{pour toute machine $i$}.}
\]
Cet algorithme entre dans le cadre \guil{map-reduce}, comme précédemment: la mise \`a jour des $(\w_i)_{k}$ et $(\u_i)_k$ correspond au map et celle de $\v_{k}$ au reduce. L'opération coûteuse qui dépend des données (le calcul des $\prox_{f_i/n\rho}$) est ainsi effectuée efficacement en parallèle.


%======================================================================
\paragraph{Bilan sur les algorithmes distribués}

La distribution et la parallélisation des calculs dans les algorithmes d'optimisation du premier ordre se font naturellement dans le cadre d'une plateforme de calcul centralisée. On peut aussi étendre ces algorithmes à des architectures plus évoluées (avec des versions asynchrones ou décentralisées); les contraintes sur les communications entre machines sont alors importantes à prendre en compte. Soulignons que ces questions sur l'adaptation des algorithmes d'optimisation à l'hétérogénéité des plateformes de calcul sont au c{\oe}ur de défis actuels en optimisation pour les big data.

%Les deux types d'algorithmes admettent différents équilibres en terme de coup par itération et de gain par itération.
%Intuitivement, les méthodes incrémentales utilisent l'information des données plus efficacement que les méthodes complètes.
%SG le coût de chaque itération est indépendant de $N$.
%mais la convergence demande $\ga_t\rightarrow 0$ 
%
%G cout est linéaire en $n$ mais pas de descente constant (ou bien choisi).


%======================================================================
%======================================================================
\subsection{Au-delà de ce chapitre}
%======================================================================
%======================================================================

L'optimisation et ses interactions avec les autres disciplines (statistique, apprentissage, informatique, calcul distribué) forment un domaine, riche et actif, de recherches et de développements, au c{\oe}ur des méthodes numériques pour l'analyse de données. 
Ce chapitre se veut un guide de lecture pour se frayer un chemin parmi les multiples documents et publications sur ce domaine. Il présente des éléments sur la formulation des problèmes d'optimisation en apprentissage et sur trois aspects des méthodes des traitement de grands volumes de données (algorithmes du premier ordre, incrémentaux, distribués). Le choix a été fait d'insister sur les idées générales, au détriment des détails (mathématiques ou algorithmiques), d'exemples concrets et des spécificités de certaines applications. 

%Pour aller plus loin que cette courte introduction, voici pour finir quelques références et des liens vers des exercices.

%\begin{itemize}
%\item \textit{\textbf{Sur les épaules des géants}}
%Pour approfondir les notions présentées dans ce chapitre, on peut se référer
%% aux nombreux documents, notamment des cours, disponibles sur internet, ainsi qu'
%aux livres ou articles de revue suivants, traitant notamment les sujets associés: [Boyd] modélisation en optimisation, [Bottou]  optimisation pour l'apprentissage, [Hull mou] analyse convexe, [Nesterov] algorithme du premier ordre, [Bubeck] algorithme pour l'apprentissage, [Wotao] optimisation distribuée.
%
%
%\item \textit{\textbf{C'est en forgeant qu'on devient forgeron}}
%Pour manipuler les notions présentées dans ce chapitre, des exercices variés sont indispensables, mais les réunir ici ne semblent pas adapté à la diversité des lecteurs et aux multiples aspects (théorique, algorithtmiques, applicatifs...) qu'il faudrait couvrir. On peut trouver des exercices de tout type (modélisation, maths, programmation, utilisation de logiciel) sur la page internet 
%
%\texttt{ljk.imag.fr/membres/Jerome.Malick/teaching.html}
%
%Ces exercices ont été testés par les étudiants de l'école d'ingénieurs ENSIMAG ou les étudiants du parcours "science des données" commun aux masters de mathématiques appliquées et d'informatique de Université Grenoble Alpes. On trouve en particulier des exercices de programmation sur la régression logistique régularisée et sur les systèmes de recommendation.
%\end{itemize}

%======================================================================
%\paragraph{Sur les épaules des géants} 
Pour approfondir les notions présentées dans cette courte introduction, on peut se référer aux livres ou articles suivants : %, traitant notamment les sujets associés: 
\citep{boyd04} pour la modélisation en optimisation, \citep{bottou17} pour une vue d'ensemble de l'optimisation pour l'apprentissage, \citep{hiriart01} pour l'analyse analyse convexe, \citep{nesterov13} pour les algorithmes du premier ordre, \citep{bubeck15} pour l'étude des algorithmes d'optimisation pour l'apprentissage.  % [nous ?] pour l'optimisation distribuée.


%======================================================================
%======================================================================
\section{Exercices}
%======================================================================
%======================================================================

Pour manipuler les notions présentées dans ce chapitre, des exercices variés sont indispensables, mais les réunir ici ne semble pas adapté à la diversité des lecteurs et aux multiples aspects (théoriques, algorithmiques, applicatifs) qu'il faudrait couvrir. Nous donnons donc uniquement quelques exercices qui précisent ou illustrent des points techniques présentés dans ce chapitre introductif. Nous renvoyons aux ouvrages cités précédemment, qui apportent des exercices et des compléments adaptés aux attentes spécifiques des lecteurs de cet vue d'ensemble du domaine.



\begin{itemize}[itemindent=20mm]
%=======
\Quest{5.1}\textbf{Opérations qui préservent la convexité.} Montrer dans les trois cas suivants, que la fonction $f$ est convexe:
\begin{itemize}
\item[a)] une somme positive de fonctions convexes $f_1$ et $f_2$
\[
f (\w) = \al_1f_1(\w) + \al_2f_2(\w) \quad \text{avec $\al_1,\al_2 \geq 0$}
\]
\item[b)] une composition d'une fonction convexe $g$ avec une application affine
\[
f(\w) = g(X\w + b)
\]
\item[c)] un maximum d'une famille de fonctions convexes $(f_i)_i$
\[
f(\w) = \max_i{f_i (\w)}
\]
\end{itemize}

%=======
\Quest{5.2}\textbf{Exemple simple de l'algorithme de gradient.} 
Soit la fonction simple $f(\om)=\om^2$ d'une variable réelle $\om\in\RR$. \'Ecrire l'itération de l'algorithme du gradient à pas constants avec cette fonction. Regarder les pas $\ga_k=1$, $\ga_k=2$ et $\ga_k=1/2$; pour lequel(s) a-t-on convergence ?


%=======
\Quest{5.3}\textbf{Opérateur proximal pour la norme $\norm{\cdot}_1$.} Nous allons montrer\index{opérateur proximal} que la fonction valeur absolue, non dérivable en $0$, est prox-simple, puis en déduire que la norme $\norm{\cdot}_1$ l'est aussi. L'opérateur proximal associé  est appelé \guil{soft-thresholding}.\index{opérateur proximal!soft-thresholding}
\begin{itemize}
\item[a)]
En explicitant les trois cas ($z=0$, $z>0$ et $z<0$), montrer que l'unique solution de \eqref{eq:prox} pour la valeur absolue s'écrit: 
\[
\prox_{|\cdot|} (\om)= \accol{ \begin{array}{cl}0 &\text{si $-1\leq\om\leq1$}\\ \om-1 &\text{si $\om>1$}\\ \om+1 &\text{si $\om<-1$} \end{array}}
\]
\item[b)] En déduire que la norme $\norm{\cdot}_1$ est prox-simple aussi, avec :
\begin{equation}\label{eq:soft}
\prox_{\norm{\cdot}_1} (\w) = (\prox_{|\cdot|} (\w_i))_{i=1,\ldots,d} ~\in \RR^d
\end{equation}
On voit ainsi que cet opérateur a pour effet de produire des vecteurs avec des composantes nulles (en annulant $\w_i\in [-1,1]$). Cet effet explique l'intérêt de la régularisation par la norme 1 pour sélectionner les attributs les plus importants.
\end{itemize}

%=======
\Quest{5.4}\textbf{\'Eclatement de prox.}  
Cet exercice propose de calculer les opérateurs proximaux\index{opérateur proximal} des deux fonctions utilisées à la section \ref{sec:distrib} pour décomposer les calculs selon chaque bloc de données. Rappelons que le point proximal est l'unique solution du problème \eqref{eq:prox}, ce qui, pour une fonction $h\colon(\RR^d)^n\rightarrow\RR$
et un point $(\w_1,\ldots,\w_n)\in (\RR^d)^n$, s'écrit :
\[%begin{equation}\label{eq:prox}
\min_{(\z_1,\ldots,\z_n)\in (\RR^d)^n} ~~h\big((\z_1,\ldots,\z_n)\big) +\frac{1}{2}\sum^n_{i=1}\norm{\z_i-\w_i}^2
\]%end{equation}


\begin{itemize}
\item[a)] Pour la fonction $f$ de \eqref{eq:feclat}, montrer que :
\[
\prox_{f}\big((\w_1,\ldots,\w_n)\big) = \Big(\prox_{f_i/n}(\w_i)\Big)_{i=1,\ldots,n}
\]
\item[b)] Pour la fonction $g$ de \eqref{eq:geclat} et en supposant que $r$ est différentiable, montrer que les composantes 
de $\prox_{g}\big((\w_1,\ldots,\w_n)\big)$ sont toutes identiques et égales à :
\[
\prox_{r/n}\left(\frac{1}{n}\sum^n_{i=1}\w_i\right)
\]
Indication: avec la différentiabilité de $r$, écrire les conditions d'optimalité \eqref{eq:condopt} (nécessaires et suffisantes, car $r$ est convexe). Le cas non différentiable se traite de la même manière avec l'aide du sous-différentiel \citep{hiriart01}.
\end{itemize}

%=======
\Quest{5.5}\textbf{Exercices sur machine.} On trouve des exercices de tout type (modélisation, maths, programmation, utilisation de logiciels...) sur la page internet \texttt{ljk.imag.fr/membres/Jerome.Malick/teaching.html}. Ils ont été testés par les étudiants de l'école d'ingénieurs ENSIMAG et ceux du parcours \guil{science des données} commun aux masters de mathématiques appliquées et d'informatique de Université Grenoble Alpes. 
En particulier, les lecteurs de ce chapitre pourraient être intéressés par les \guil{jupyter notebooks} avec des exercices sur la régression logistique régularisée et sur les systèmes de recommendation.
\end{itemize}





%%======================================================================
%\bibliographystyle{plainnat-fr}
%\bibliography{biblio-Optim}


%%======================================================================
%\paragraph{} 



%%======================================================================
%%======================================================================
%\section{Other algos !}
%%======================================================================
%%======================================================================
%
%%======================================================================
%\paragraph{traiter une variable à la fois}
%
%\[
%\min_{\w} f(\w)=f(\w_1,\ldots,\w_N)
%\]
%la variable $\w\in \RR^d$ est partionnée en $N$ blocs de variables.
%
%les algorithmes d'optimisation par (blocs de) coordonées consistent à mettre à jour les blocs de variable un par un en gardant les autres fixes. Dans le cas d'une etape de gradient partiel
%\[
%x^{k+1}_i = x^{k}_i - \ga_t \nabla_i f(\w^k) \quad \text{et} \quad x^{k+1}_j = x^{k}_j \text{pour $j\neq i$}
%\]
%
%$\Om$ est convexe non-lisse séparable
%
%\[
%\min_\w f(x_1,\ldots,x_n) + \la\sum^m_{i=1} \Om_i(x_i)
%\]
%mise a jour des variables $x_i$ une par une, de manière cyclique, aléatoire, ou parallèle. Dans de nombreux cas, 
%ces méthodes s'avèrent être plus rapide que les méthodes qui mettent à jour tous les variables. Insistons sur le fait que cette approche est limité au cas où la partie non-lisse est séparable.
%
%
%%======================================================================
%\paragraph{algorithmes duaux}
%
%problème primal : traiter un bloc de données à la fois correspond à traiter une seule fonction $f_i$ à la fois.
%
%\[
%f_i(\trans{A_i}\w)
%\]
%
%problème dual : traiter un bloc de données à la fois correspond à traiter une variable $u_i\in \RR^m$ à la fois.
%
%\[
%\max_{u = (u_1,\ldots,u_n)\in \RR^{m\cdot n}} \theta(u) = 
%-\la \Om^*\Big(\frac{1}{\la n}\sum^n_{i=1} A_iu_i \Big) - \frac{1}{n}\sum^n_{i=1}f_i^*(-u_i)
%\]
%
%Un exemple très connu en apprentissage est la classification par \guill{SVM} (support vector machines, en anglais) où le problème primal est
%\[
%\min_\w \sum^n_{i=1} \max\{0,1-y_i\trans{a_i}\w\} + \frac{1}{2}{\norm{\w}_2}^2
%\]
%C'est un probleme fortement convexe par la présence de la régularisation par norme $\ell_2$; son dual (lisse) prend la forme d'un problème quadratique avec contrainte de boite sur la variable.
%\[
%\min_{0\leq u \leq 1/\la} \frac{1}{2}\trans{u}A\trans{A}u - \sum^n_{i=1} u_i
%\]

%%======================================================================
%%======================================================================
%\section{Recomm}
%%======================================================================
%%======================================================================
%
%%======================================================================
%\paragraph{Exemple: systèmes de recommendation}
%
%
%dont l'exemple le plus connu est la recommendation de films à un utilisateur de service de film en flux continu via Internet
%
%les données paire (utilisateur, film)
%
%Une modélisation de ce problème est de fixer le rang $k$ de la matrice $X$; on cherche alors les deux facteurs 
%$P \in \RR^{m\times k},Q\in \RR^{k\times n}$ tel que $X=\trans{P}Q$ nous permet d'apprendre au mieux les données. Le modèle d'apprentissage est complété par une régularisation par norm $\ell_2$ sur les matrices. Cette régularisation est particulièrement intéressante dans le cas car, en plus d'apporter un meilleur conditionneent 
% 
%la factorisation de matrice 
%\[
%\min_{P,Q} \sum_{i,j: \text{$r_{ij}$ exists}}(r_{ij}-\ip{p_i}{q_j})^2 +\frac{\la}{2}(\sum_{i}\norm{p_i}^2 + \frac{\la}{2}(\sum_{j}\norm{q_j}^2 
%\]
%
%
%
%
%
%Exemple : gradient incrémental pour les systèmes de recommendation
%
%On choisit aléatoirement un couple (utilisateur, film) et 
%on calcule le gradient partiel. Il se trouve que le gradient ne concerne que les couple $(i,j)$ -- ce qui est particulier à cet exemple et n'est pas vrai en général. Des calculs élémentaires mais entravé d'indices donnent en effet l'expression
%\[
%\nabla f_{i,j}(p, q) = \left[\begin{array}{c}
%-2(y_{ij} - \trans{p_i} q_{j}) q_{j} + \la p_i\\
%-2(y_{ij} - \trans{p_i} q_{j}) p_{i} + \la q_j
%\end{array}\right]
%\]
%Dans ce cas particulier, une itération du gradient stochastique s'écrit ainsi
%\[
%\left[\begin{array}{c}
%p_{i,k+1}\\
%q_{i,k+1}
%\end{array}\right] = 
%\left[\begin{array}{c}
%p_{i,k}\\
%q_{i,k}
%\end{array}\right]  -\ga \nabla f_{i,j}(p_k, q_k)
%\]
%



%
%
%%======================================================================
%\section{Exercices}
%
%exo d'accompagnemen du cours sur les notions de base.
%
%des exos specifiques sur des applications en analyse de données.
%
%idealement mélangés comme en 1
%
%\paragraph{exercice}
%
%comment reconnaitre la convexité ?
%
%
%Montrer la convexité de fonction des moindres carrés
%Montrer la convexité de la fonction des SVM.
%
%
%exo sur les fonctions diff ?
%
%
%exo sur le gradient a pas constant ?
%
