Eclats de vers : Matemat 10 : Optimisation - 2

Soit un espace vectoriel \(E\) sur \(\setR\) et des vecteurs \(x,y \in E\), où \(x \ne 0\). La projection de \(y\) sur \(x\) est le vecteur :

\[p \in \combilin{x} = \{ \gamma \cdot x : \gamma \in \setR \}\]

qui minimise sur \(\combilin{x}\) la norme usuelle \(\norme{e} = \sqrt{ \scalaire{e}{e} }\) de l'écart \(e\) séparant \(y\) de \(p\). Afin de trouver le \(\lambda \in \setR\) minimisant cet écart, on utilise la fonction \(\mathcal{E} : \setR \mapsto \setR\) définie par :

\[\mathcal{E}(\gamma) = \scalaire{y - \gamma \cdot x}{y - \gamma \cdot x} = \scalaire{y}{y} - 2 \cdot \gamma \cdot \scalaire{x}{y} + \gamma^2 \cdot \scalaire{x}{x}\]

pour tout \(\gamma \in \setR\). Nous imposons l'annulation de la dérivée en un certain \(\gamma = \lambda\) qui minimise potentiellement l'écart :

\[\partial \mathcal{E}(\lambda) = - 2 \cdot \scalaire{x}{y} + 2 \cdot \lambda \cdot \scalaire{x}{x} = 0\]

ce qui nous donne :

\[\lambda = \frac{\scalaire{x}{y}}{\scalaire{x}{x}}\]

Calculons à présent la norme de \(y\) en partant de la relation :

\[y = p + e\]

On déduit de la propriété d'orthogonalité que :

\[\scalaire{p}{e} = \lambda \cdot \scalaire{x}{e} = 0\]

On peut donc appliquer le théorème de pythagore, qui nous dit que :

\[\norme{y}^2 = \norme{p + e}^2 = \norme{p}^2 + \norme{e}^2\]

Si :

\[D = \min_{\gamma \in \setC} \norme{y - \gamma \cdot x}\]

on a donc :

\[D^2 = \norme{y - \lambda \cdot x}^2 = \norme{e}^2 = \norme{y}^2 - \norme{p}^2\]

et finalement :

\[D^2 = \norme{y}^2 - \abs{\lambda}^2 \cdot \norme{x}^2\]

Par positivité du produit scalaire, on a bien évidemment :

\[\norme{e}^2 = \scalaire{e}{e} \ge 0\]

Le théorème de Pythagore implique donc que :

\[\scalaire{p}{p} = \norme{p}^2 \le \norme{y}^2 = \scalaire{y}{y}\]

En substituant \(p = \lambda \cdot x\), on obtient :

\[\conjaccent \lambda \cdot \lambda \cdot \scalaire{x}{x} \le \scalaire{y}{y}\]

Mais comme :

\[\conjaccent{\lambda} = \conjugue \frac{ \scalaire{x}{y} }{ \scalaire{x}{x} } = \frac{ \scalaire{y}{x} }{ \scalaire{x}{x} }\]

on a :

\[\frac{ \scalaire{y}{x} }{ \scalaire{x}{x} } \cdot \frac{ \scalaire{x}{y} }{ \scalaire{x}{x} } \cdot \scalaire{x}{x} \le \scalaire{y}{y}\]

En simplifiant et en multipliant par la norme carrée de \(x\), on a finalement :

\[\scalaire{y}{x} \cdot \scalaire{x}{y} \le \scalaire{x}{x} \cdot \scalaire{y}{y}\]

relation dont la racine nous donne l'inégalité de Cauchy-Schwartz :

\[\abs{ \scalaire{x}{y} } \le \norme{x} \cdot \norme{y}\]

L'égalité :

\[\norme{p}^2 = \norme{y}^2\]

n'est atteinte que lorsque \(\norme{e}^2 = 0\), c'est-à-dire \(e = 0\) et \(y = \lambda \cdot x\) pour un certain \(\lambda \in \setC\). Ce constat nous amène aux problèmes d'optimisations suivants. Soit un vecteur \(c \ne 0\) fixé et :

\[d = \norme{c} = \sqrt{ \scalaire{c}{c} }\]

On cherche à maximiser ou à minimiser :

\[\varphi(x) = \scalaire{c}{x}\]

sur l'ensemble \(B = \boule(0,r) = \{ x : \norme{x} \le r \}\). L'inégalité de Cauchy-Schwartz nous dit que :

\[- d \cdot r \le - \norme{c} \cdot \norme{x} \le \scalaire{c}{x} \le \norme{c} \cdot \norme{x} \le d \cdot r\]

Nous allons chercher nos solutions sous la forme \(x = \alpha \cdot c\), pour un certain \(\alpha \in \setC\). On a alors :

\[\norme{\alpha \cdot c} = \abs{\alpha} \cdot \norme{c} = \abs{\alpha} \cdot d \le r\]

et donc \(\abs{\alpha} \le r / d\). Si on choisit \(\alpha = r / d\), on a :

\[\scalaire{c}{x} = \frac{r}{d} \cdot \scalaire{c}{c} = \frac{r}{d} \cdot d^2 = d \cdot r\]

La borne supérieure est atteinte. La fonction \(\varphi\) est donc maximisée :

\[\eta = \frac{r}{d} \cdot c \in \arg\max_{x \in B} \scalaire{c}{x}\]

Si on choisit \(\alpha = - r / d\), on a :

\[\scalaire{c}{x} = - \frac{r}{d} \cdot \scalaire{c}{c} = - d \cdot r\]

La borne inférieure est atteinte. La fonction \(\varphi\) est donc minimisée :

\[\theta = - \frac{r}{d} \cdot c \in \arg\min_{x \in B} \scalaire{c}{x}\]

Soit l'espace vectoriel \(E\) sur \(\setR\) et un sous-espace vectoriel \(U \subseteq E\) possédant une base orthonormée \((u_1,...,u_n)\). La projection d'un vecteur quelconque \(z \in E\) sur \(U\) est le vecteur \(p \in U\) qui minimise la norme de l'écart entre \(z\) et \(p\). On cherche donc le \(\lambda = (\lambda_1,...,\lambda_n) \in \setR^n\) qui minimise la fonction \(\mathcal{E} : \setR^n \mapsto \setR\) définie par :

\[\mathcal{E}(\gamma) = \norme{z - \sum_{i = 1}^n \gamma_i \cdot u_i}^2\]

pour tout \(\gamma = (\gamma_1,...,\gamma_n) \in \setR^n\). En utilisant \(\scalaire{u_i}{u_j} = \indicatrice_{ij}\), on obtient :

Imposant \(\partial_k \mathcal{E}(\lambda_1,...,\lambda_n) = 0\), on obtient les \(n\) équations :

\[-2 \lambda_k \cdot \scalaire{u_k}{z} + 2 \lambda_k = 0\]

On en déduit que le choix :

\[\lambda = (\lambda_1,...,\lambda_n) = (\scalaire{u_1}{z},...,\scalaire{u_n}{z})\]

minimise potentiellement \(\mathcal{E}\). Notre projection potentielle de \(z\) sur \(U\) est donc donnée par :

\[p = \sum_i \scalaire{u_i}{z} \cdot u_i\]

On peut réécrire la projection de \(z\) sur \(U\) sous la forme :

\[p = \sum_i u_i \cdot \scalaire{u_i}{z}\]

Cette expression ressemble à la contraction d'ordre \(1\) d'un tenseur d'ordre \(2\). Effectivement, si on pose :

\[\mathcal{P} = \sum_i u_i \otimes u_i\]

on a alors :

\[p = \mathcal{P} \cdot z = \contraction{\mathcal{P} }{1}{z}\]

Nous allons construire une suite de vecteurs orthonormaux \((u_1,u_2,...u_n)\) à partir d'une suite \((a_1,a_2,...,a_n)\) de vecteurs linéairement indépendants de \(E\). L'indépendance linéaire nous garantit que \(a_1 \ne 0\). On peut donc normaliser pour obtenir le premier vecteur de la série :

\[u_1 = \frac{a_1}{ \norme{a_1} }\]

On a alors :

\[\scalaire{u_1}{u_1} = \frac{ \scalaire{a_1}{a_1} }{ \norme{a_1}^2 } = \frac{ \scalaire{a_1}{a_1} }{ \scalaire{a_1}{a_1} } = 1\]

Nous projetons \(a_2\) sur \(u_1\) en utilisant le tenseur \(\mathcal{P}_1 = u_1 \otimes u_1\) et nous évaluons l'écart :

\[e_2 = (\tenseuridentite - \mathcal{P}_1) \cdot a_2 = a_2 - u_1 \cdot \scalaire{u_1}{a_2}\]

Les propriétés d'orthogonalité de l'écart nous assurent alors que \(\scalaire{u_1}{e_2} = 0\). L'indépendance linéaire nous garantit que \(e_2 \ne 0\). On peut donc normaliser en utilisant \(u_2 = e_2 / \norme{e_2}\). On a alors clairement \(\scalaire{u_2}{u_2} = 1\) et \(\scalaire{u_1}{u_2} = 0\).

Supposons à présent avoir trouvé la suite orthonormée \((u_1,u_2,...,u_k)\), où \(k \le n - 1\). Nous projetons \(a_{k + 1}\) sur l'espace vectoriel \(U_k = \combilin{u_1,u_2,...,u_k}\) en utilisant le tenseur :

\[\mathcal{P}_k = \sum_{i = 1}^k u_i \otimes u_i\]

et nous évaluons l'écart :

\[e_{k + 1} = (\tenseuridentite - \mathcal{P}_k) \cdot a_{k + 1} = a_{k + 1} - \sum_{i = 1}^k u_i \cdot \scalaire{u_i}{a_{k + 1} }\]

Les propriétés d'orthogonalité de l'écart nous assurent alors que \(\scalaire{u_j}{e_{k + 1} } = 0\) pour tout \(j \in \{1,2,...,k\}\). L'indépendance linéaire nous garantit que \(e_{k + 1} \ne 0\). On peut donc normaliser en utilisant :

\[u_{k + 1} = \frac{e_{k + 1} }{ \norme{ e_{k + 1} } }\]

On a alors clairement \(\scalaire{u_j }{u_{k + 1} } = \indicatrice_{j, k + 1}\).

Nous disposons donc à la fin du processus d'une suite de vecteurs \((u_1,u_2,...,u_n)\) orthonormée :

\[\scalaire{u_i}{u_j} = \indicatrice_{ij}\]

Ce procédé est appelé processus d'orthogonalisation de Gram-Schmidt.

Soient \(u_1,u_2,...,u_m \in \matrice(\corps,n,1)\) des vecteurs matriciels orthonormés pour le produit scalaire usuel :

\[u_i^\dual \cdot u_j = \indicatrice_{ij}\]

La matrice de projection associée au tenseur de projection sur \(U = \combilin{u_1,...,u_m}\) s'écrit :

\[P = \sum_{k = 1}^m u_k \otimes u_k = \sum_{k = 1}^m u_k \cdot u_k^\dual\]

Il s'agit donc d'une matrice de taille \((n,n)\). Elle permet de projeter tout vecteur matriciel \(z \in \matrice(\corps,n,1)\) sur \(U\) :

\[P \cdot z = \sum_{k = 1}^m u_k \cdot u_k^\dual \cdot z = \sum_{k = 1}^m u_k \cdot \scalaire{u_k}{z}\]

En terme de composantes, si \(u_k = ( u_{kj} )_j\) et si \(P = ( p_{ij} )_{i,j}\), on a :

\[p_{ij} = \sum_{k = 1}^m u_{ki} \cdot \conjaccent{u}_{kj}\]

expression qui ressemble furieusement à un produit matriciel. Soit la matrice \(U\) de taille \((n,m)\) rassemblant les \(m\) vecteurs \(u_k\) :

\[U = [u_1 \ u_2 \ \hdots \ u_m]\]

On a alors :

Le produit ci-dessus peut donc s'écrire :

\[p_{ij} = \sum_{k = 1}^m b_{ik} \cdot c_{kj}\]

où \(b_{ik} = \composante_{ik} U\) et \(c_{kj} = \composante_{kj} U^\dual\). On a donc finalement :

\[P = U \cdot U^\dual\]

Soit une matrice de produit scalaire \(A \in \matrice(\corps,n,n)\). On peut utiliser le procédé de Gram-Schmidt pour construire une base de vecteurs orthogonaux pour le produit scalaire :

\[\scalaire{x}{y} = x^\dual \cdot A \cdot y\]

On part d'une suite \((a_1,a_2,...,a_n)\) de vecteurs matriciels linéairement indépendants de \(\corps^n\) (typiquement la base canonique). On commence par normaliser \(a_1\) :

\[u_1 = \frac{a_1}{ \sqrt{a_1^\dual \cdot A \cdot a_1} }\]

et on construit étape après étape la suite des \(u_i\). Supposons être arrivé à la suite \((u_1,u_2,...,u_k)\) vérifiant \(\scalaire{u_i}{u_j} = u_i^\dual \cdot A \cdot u_j = \indicatrice_{ij}\), où \(k \le n - 1\). Nous projetons \(a_{k + 1}\) sur l'espace vectoriel \(\combilin{u_1,u_2,...,u_k}\) et nous évaluons l'écart :

Ensuite, nous normalisons le résultat :

\[u_{k + 1} = \frac{e_{k + 1} }{ \sqrt{e_{k + 1}^\dual \cdot A \cdot e_{k + 1} } }\]

Nous disposons donc à la fin du processus d'une suite de vecteurs \((u_1,u_2,...,u_n)\) orthonormée :

\[\scalaire{u_i}{u_j} = u_i^\dual \cdot A \cdot u_j = \indicatrice_{ij}\]

On dit que la suite des \(u_i\) est $A$-orthonormée. Si on définit la famille de matrices :

\[U_m = [u_1 \ u_2 \ ... \ u_m]\]

pour tout \(m \le n\), on peut réecrire l'orthogonalité comme suit :

\[U_m^\dual \cdot A \cdot U_m = I_m\]

On note aussi \(U = U_n\).

Nous allons présenter des algorithmes permettant de résoudre approximativement des problèmes de minimisation d'une fonction \(\varphi\) sur \(\setR^n\). Ces algorithmes partent d'un point initial \(x_0 \in \Omega\) et itèrent schématiquement comme suit :

\[x_{k + 1} = I(x_k) = x_k + p_k\]

pour un certain \(p_k \in \setR^n\). On espère bien entendu que la suite converge et que :

\[x_N \approx \lim_{k \to \infty} x_k = \arg\min_{x \in \Omega} \varphi(x)\]

pour \(N\) assez grand. Nous adoptons les notations :

\[J = (\partial \varphi)^\dual\]

pour le gradient, de taille \((n,1)\), et :

\[H = \partial^2 \varphi\]

pour la hessienne, de taille \((n,n)\). On note également :

\[\Phi_k = \Phi(x_k)\]

pour toute fonction \(\Phi\) (par exemple, \(\Phi \in \{\varphi,J,H\}\)).

Soit l'itération :

\[x_{k+1} = x_k - \alpha_k \cdot p_k\]

où \(\alpha_k \in \setR\) et \(p_k \in \setR^n\). On choisit généralement le paramètre \(\alpha_k\) de façon à minimiser le développement d'ordre deux :

\[\varphi_{k + 1} \approx \varphi_k - \alpha_k \cdot J_k^\dual \cdot p_k + \frac{\alpha_k^2}{2} \cdot p_k^\dual \cdot H_k \cdot p_k\]

En imposant l'annulation de la dérivée de ce développement par rapport à \(\alpha_k\), on en déduit que :

\[- J_k^\dual \cdot p_k + \alpha_k \cdot p_k^\dual \cdot H_k \cdot p_k = 0\]

La valeur optimale de \(\alpha_k\) s'écrit :

\[\alpha_k = \frac{J_k^\dual \cdot p_k}{p_k^\dual \cdot H_k \cdot p_k}\]

La méthode de la plus grande pente consiste à partir à chaque itération du point \(x^{(k)}\) et à suivre la direction \(\delta_k\) où \(\varphi\) descend le plus rapidement dans le voisinage immédiat. En première approximation, si :

\[x_{k + 1} = x_k + \alpha_k \cdot \delta_k\]

pour un certain \(\alpha_k \in \setR\), on a :

\[\varphi_{k+1} \approx \varphi_k + \alpha_k \cdot J_k^\dual \cdot \delta_k\]

Nous choisissons le vecteur \(\delta_k\) qui minimise \(J_k^\dual \cdot \delta_k = \scalaire{J_k}{\delta_k}\) sur \(\boule(0,\norme{J_k})\), c'est-à-dire :

\[\delta_k = - J_k\]

On a alors :

\[x_{k + 1} = x_k - \alpha_k \cdot J_k\]

La valeur optimale de \(\alpha_k\) s'écrit donc :

\[\alpha_k = \frac{J_k^\dual \cdot J_k}{J_k^\dual \cdot H_k \cdot J_k}\]

Il s'agit ici d'optimiser le pas \(s_k\) :

\[x_{k + 1} = x_k + s_k\]

pour minimiser le développement :

\[\varphi_{k+1} \approx \varphi_k + J_k^\dual \cdot s_k + \unsur{2} \cdot s_k^\dual \cdot H_k \cdot s_k\]

Mais comme \(H = H^\dual\), l'annulation du gradient par rappord à \(s_k\) nous donne :

\[J_k + H_k \cdot s_k \approx 0\]

On en déduit la valeur optimale :

\[s_k = - H_k^{-1} \cdot J_k\]

Soit l'itération :

\[x_{k + 1} = x_k - \alpha_k \cdot p_k\]

où \(\alpha_k \in \setR\) et \(p_k \in \setR^n\). On a comme d'habitude :

\[\alpha_k = \frac{J_k^\dual \cdot p_k}{p_k^\dual \cdot H_k \cdot p_k}\]

Lorsque \(k = 0\), on prend le gradient comme direction :

\[p_0 = J_0\]

Pour \(k \ge 1\), on choisit \(p_k\) comme combinaison linéaire du gradient \(J_k\) et du pas précédent \(p_{k - 1}\) :

\[p_k = J_k - \beta_k \cdot p_{k-1}\]

où \(\beta_k \in \setR\). L'idée des gradients conjugué est de construire une suite de \(p_k\) orthogonaux entre-eux afin de minimiser la fonction dans toutes les directions. Au bout de \(n\) itérations, on espère avoir construit une base de \(\setR^n\) et être très proche du minimum global. En fait, nous n'allons pas vérifier que toutes les directions sont orthogonales, mais seulement que deux directions successives le sont. On demande donc l'orthogonalité au sens du produit scalaire défini par \(H\) :

\[p_k^\dual \cdot H_k \cdot p_{k - 1} = J_k^\dual \cdot H_k \cdot p_{k - 1} - \beta_k \cdot p_{k - 1}^\dual \cdot H_k \cdot p_{k - 1} = 0\]

On en déduit la valeur de \(\beta_k\) :

\[\beta_k = \frac{J_k^\dual \cdot H_k \cdot p_{k - 1}}{p_{k - 1}^\dual \cdot H_k \cdot p_{k - 1}}\]

Nous allons à présent obtenir une valeur approximative de \(\beta_k\) en fonction des variations du gradient. Le développement d'ordre un de \(J\) s'écrit :

\[J_k - J_{k - 1} \approx H_k \cdot (x_k - x_{k - 1}) = - \alpha_k \cdot H_k \cdot p_{k - 1}\]

On en déduit que :

\[\beta_k \approx \frac{J_k^\dual \cdot (J_k - J_{k - 1})}{p_{k - 1}^\dual \cdot (J_k - J_{k - 1})}\]

Soit la fonction \(f : \setR^n \mapsto \setR^m\). On cherce à minimiser la fonction \(\varphi = f^\dual \cdot f / 2\). Le problème de minimisation s'écrit alors :

\[\arg\min_x \unsur{2} f(x)^\dual \cdot f(x)\]

Dans ce cas, si on définit :

\[D = \partial f\]

le gradient s'écrit :

\[J = D^\dual \cdot f\]

Si on suppose que les dérivées secondes de \(f\) sont négligeables par rapport aux dérivées premières, on a :

\[H \approx D^\dual \cdot D\]

La méthode de Newton devient dans ce cas :

\[x_{k+1} = x_k + s_k\]

avec :

\[s_k = -[D_k^\dual \cdot D_k]^{-1} \cdot D_k^\dual \cdot f_k\]

C'est une variante de la méthode des moindres carrés, utile dans les cas où \(D_k^\dual \cdot D_k\) est numériquement proche d'une matrice non inversible. On ajoute alors la matrice identité multipliée par un scalaire \(\lambda\) sur la diagonale :

\[s_k = -[D_k^\dual \cdot D_k+ \lambda_k \cdot I]^{-1} \cdot D_k^\dual \cdot f_k\]