A x = b vu comme un problème d'optimisation

Pour résoudre $A {\bf x} = {\bf b}$ on va chercher le minimum de la fonctionnelle

$$ J({\bf x}) = \frac{1}{2} \, {\bf x}^T \, A \, {\bf x} - {\bf b}\, . {\bf x} $$

En ce point la dérivée s'annule et dans un cas précis la dérivée est justement $A {\bf x} - {\bf b}$.

Calcul de dérivée

La notion de dérivée en dimension supérieure à 1 doit être manipulée avec précaution car il y a différents types de dérivées, les dérivées partielles et la dérivée totale.

Ce qui nous intéresse est la dérivée dans une direction ce qui correpond à la dérivée partielle en $y$ si on va dans la direction de l'axe $y$.

Définition

On dit que $f : \Omega \subset X \rightarrow Y$ ($\Omega$ ouvert) est dérivable en ${\bf p} \in \Omega$ si

$$ \exists \; f'({\bf p}) \in L(X,Y) \textrm{ tel que} \\ f ({\bf p} + {\bf h}) = f({\bf p}) + f'({\bf p}) ( {\bf h} ) + {\bf h} \; \varepsilon({\bf h}) \quad \textrm{ avec } \lim_{\bf h \rightarrow 0} \varepsilon({\bf h}) = {\bf 0} $$

où $L(X,Y)$ représente l'espace des applications linéraires continues de $X$ dans $Y$. Attention c'est $f'(p) \in L(X,Y)$ et non $f'$.

On note que si $f$ est dérivable en ${\bf p}$ alors $\forall \; {\bf h} \in X$

$$ f'({\bf p}) ({\bf h})= \lim_{\theta \rightarrow 0} \frac{f({\bf p} + \theta \, {\bf h}) - f({\bf p})}{\theta} $$

Attention à bien regarder le type de chaque terme. Pour bien visualiser j'ai considéré que $f$ est une fonction scalaire donc $Y = ℝ$. Par contre $X = ℝ^n$ avec $n > 1$.

Autre notation

Je trouve la notation avec le gradient plus claire. On voit mieux où sont les vecteurs et les scalaires :

$$ f ({\bf p} + {\bf h}) = f({\bf p}) + \nabla \! f\,({\bf p})^T \; {\bf h} + {\bf h}^T \; \varepsilon({\bf h}) \quad $$

Comme j'ai pris $Y = ℝ$ l'application $f$ est scalaire. $\nabla \! f \, ({\bf p})$ est un vecteur dont le produit scalaire avec ${\bf h}$ donne un réel.

Calculons la dérivée de J suivant une direction

On calcule la dérivée de $J({\bf x})$ au point ${\bf p}$ suivant la direction donnée par le vecteur $\bf h$ :

$$ \begin{align} J'({\bf p}) ({\bf h}) &= \lim_{\theta \rightarrow 0} \; \frac{J({\bf p} + \theta \, {\bf h}) - J({\bf p})}{\theta} \\ &= \lim_{\theta \rightarrow 0} \; \frac{1}{\theta} \; \left( \frac{1}{2} \; ({\bf p} + \theta {\bf h})^T \, A \; ({\bf p} + \theta {\bf h}) - {\bf b}^T \, ({\bf p} + \theta {\bf h}) - \frac{1}{2} \; {\bf p}^T \, A \; {\bf p} + {\bf b}^T \, {\bf p} \right) \\ &= \lim_{\theta \rightarrow 0} \; \frac{1}{\theta} \; \left( \frac{1}{2} \left( \theta \, {\bf p}^T \, A \; {\bf h} + \theta \, {\bf h}^T \, A \; {\bf p} + \theta^2 \, {\bf h}^T \, A \; {\bf h} \right) - \theta \, {\bf b}^T \, {\bf h} \right) \\ &= \frac{1}{2} \left( {\bf p}^T \, A \; {\bf h} +{\bf h}^T \, A \; {\bf p} \right) - {\bf b}^T \, {\bf h} \\ &= \frac{1}{2} \left( {\bf p}^T \, A \; {\bf h} +{\bf p}^T \, A^T \; {\bf h} \right) - {\bf b}^T \, {\bf h} \\ &= \frac{1}{2} \left( {\bf p}^T \, ( A + A^T ) - {\bf b}^T \right) \, {\bf h} \\ &= \frac{1}{2} \left( ( A + A^T ) \, {\bf p} - {\bf b} \right) \, . \, {\bf h} \end{align} $$

donc

$$ \begin{align} J' : \; & {\bf x} \in \Omega \subset ℝ^n \longrightarrow \;L(ℝ^n, ℝ) \\ & {\bf x} \; \longmapsto \; \frac{1}{2} ( A + A^T ) \, {\bf x}^T - {\bf b} \end{align} $$

A symétrique

Dans le cas où A est symétrique on a

$$ J'({\bf x}) = A \; {\bf x} - {\bf b} $$

On note que si la dérivée s'annule alors on a résolu notre système matriciel !

Il faut donc que

• A soit symétrique
• J ait un minimum (par certain, on a vu un cas de point selle)

et dans ce cas on peut utiliser une méthode de gradient pour trouver le minimum qui sera la solution de notre système matriciel.

Propriété

Si A est symétrique et définie positive alors $J$ est convexe strictement et coervice (c.a.d. $\lim_{||{\bf p}|| \rightarrow \infty} J({\bf p}) = + \infty$) ce qui garantit qu'elle a un minimum.

Gradient et dérivée

Le gradient est composé des dérivées partielles suivant chaque direction de la base. On peut rérouler les calculs et calculer $\nabla J$ facilement (mais c'est long). Faisons le pour une matrice sur la partie compliquée :

$$ \begin{align} \frac{1}{2} {\bf x}^T A {\bf x} &= \frac{1}{2} \begin{pmatrix} x_1 & x_2 \end{pmatrix} \begin{pmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} \\ &= \frac{1}{2} \left( x_1 a_{11} x_1 + x_1 a_{12} x_2 + x_2 a_{21} x_1 + x_2 a_{22} x_2 \right) \\ & = \frac{1}{2} \left( a_{11} x_1^2 + (a_{12} + a_{21}) x_1 x_2 + a_{22} x_2^2 \right) \end{align} $$

Calcul des dérivées partielles de $ \frac{1}{2} x^T A x $ :

$$ \frac{\partial}{\partial x_1} \left( \frac{1}{2} (a_{11} x_1^2 + (a_{12} + a_{21}) x_1 x_2 + a_{22} x_2^2) \right) = \frac{1}{2} (2 a_{11} x_1 + (a_{12} + a_{21}) x_2) = a_{11} x_1 + \frac{a_{12} + a_{21}}{2} x_2 $$

$$ \frac{\partial}{\partial x_2} \left( \frac{1}{2} (a_{11} x_1^2 + (a_{12} + a_{21}) x_1 x_2 + a_{22} x_2^2) \right) = \frac{1}{2} ((a_{12} + a_{21}) x_1 + 2 a_{22} x_2) = \frac{a_{12} + a_{21}}{2} x_1 + a_{22} x_2 $$

Donc le gradient de $\frac{1}{2} {\bf x}^T A {\bf x}$ est :

$$ \nabla \left( \frac{1}{2} x^T A x \right) = \begin{pmatrix} a_{11} x_1 + \frac{a_{12} + a_{21}}{2} x_2 \\ \frac{a_{12} + a_{21}}{2} x_1 + a_{22} x_2 \end{pmatrix} = \frac{1}{2} (A + A^T) \, {\bf x} $$

Il faut faire le calcul pour une matrice $n \times n$ mais on arrive au même résultat et donc

$$ J'({\bf x}) = \nabla J({\bf x}) $$

Ainsi on peut aussi calculer la dérivée suivant la direction $\bf h$ avec le gradient :

$$ J'({\bf x})({\bf h}) = \nabla J({\bf x}) \, . \, {\bf h} $$

import numpy as np
import plotly.offline as py
import plotly.graph_objects as go

%matplotlib inline
%config InlineBackend.figure_format = 'retina'

A = np.array([[2,1], [1,2]])
b = np.array([-1, 0.5])

# lets make the mesh
x = np.linspace(-20,20,100)
y = np.linspace(-20,20,100)
mx, my = np.meshgrid(x,y) 
M = np.array(np.meshgrid(x,y)).swapaxes(0,2)

# and compute 0.5 x.T A x - b x on each node of the mesh
MA = np.einsum("ijk, ka -> ija", M, A)
mz = 0.5 * np.einsum("ijk, ijk -> ij", MA, M) - M @ b

fig = go.Figure(data=[go.Surface(x=mx, y=my, z=mz)])
py.iplot(fig, show_link=False)

La contrainte de la propriété peut sembler être une limite trop dure mais il y a de nombreux problèmes qui génèrent des matrices A qui respectent ces contraintes. Par exemples les problèmes de simuation numérique comme celui de l'avion qu'on a vu génère une telle matrice, qui de plus sera très creuse avec une dizaine de valeur non nulle par ligne même pour une matrice de 1 milliard par 1 milliard.

C'est pour ce type de matrice qu'un méthode itérative comme le gradient offre un gain de performance incroyable puisqu'on a la solution en quasiment $O(n^2)$ au lieu de $O(n^3)$.