Dans le cursus ECG, l’algèbre linéaire en deuxième année aborde très largement le calcul des valeurs propres. Parmi ces dernières, une possède une particularité que toutes les autres n’ont pas : la valeur propre dominante. Son calcul en Python repose sur un algorithme particulier qui permet d’aboutir à une approximation convergente de la valeur de cette valeur propre toute particulière. C’est un sujet « pont », capable de tomber aussi bien en exercice d’informatique de Maths I qu’en dernière partie d’un problème de Maths II pour illustrer la convergence de suites de matrices (matrices probabilistes, par exemple). Abordons donc d’abord les fondements théoriques avant de présenter mathématiquement l’algorithme, pour voir ensuite son application en Python. Nous conclurons cet article par quelques points de vigilance (non exhaustifs malheureusement) autour de ce script.
Fondements théoriques
Soit \(A \in \mathcal{M}_n(\mathbb{R})\) une matrice carrée d’ordre \(n\). Supposons que \(A\) soit diagonalisable. Il existe donc une base de vecteurs propres \((u_1, u_2, \dots, u_n)\) associés aux valeurs propres \((\lambda_1, \lambda_2, \dots, \lambda_n)\).
Jusque-là, rien de plus classique à l’égard d’un cours de deuxième année d’algèbre !
La condition de dominance
On parle de valeur propre dominante \(\lambda_1\) lorsque sa valeur absolue est strictement supérieure à celle de toutes les autres :
\[ |\lambda_1| > |\lambda_2| \geq |\lambda_3| \geq \dots \geq |\lambda_n| \]
Cette condition est cruciale. Si deux valeurs propres ont la même valeur absolue (par exemple \(5\) et \(-5\)), l’algorithme ne convergera pas vers un point unique, mais oscillera.
La décomposition spectrale
Prenons un vecteur \(v_0\) quelconque.
Puisque \((u_i)\) est une base, il existe des scalaires \((\alpha_i)\) tels que :
\[ v_0 = \sum_{i=1}^{n} \alpha_i u_i \]
En appliquant \(k\) fois la matrice \(A\), on obtient par linéarité :
\[ A^k v_0 = \lambda_1^k \left( \alpha_1 u_1 + \sum_{i=2}^{n} \alpha_i \left(\frac{\lambda_i}{\lambda_1}\right)^k u_i \right) \]
Puisque \(|\lambda_1| > |\lambda_i|\) pour \(i > 1\), alors \(\left( \frac{\lambda_i}{\lambda_1} \right)^k \underset{k \to +\infty}{\longrightarrow} 0\). La direction de \(A^k v_0\) s’aligne donc sur celle du vecteur propre \(u_1\).
L’algorithme de la puissance
L’objectif de l’algorithme est d’extraire numériquement le couple \((\lambda_1, u_1)\) grâce au procédé qui suit.
La normalisation
Si l’on calcule simplement \(v_{k+1} = A v_k\), les composantes peuvent exploser (si \(|\lambda_1| > 1\)) ou s’annuler (si \(|\lambda_1| < 1\)). La solution consiste à normaliser le vecteur (c’est-à-dire à le diviser par sa norme) à chaque étape (nous retrouverons ce procédé dans le script ci-dessous) :
\[ v_{k+1} = \frac{A v_k}{\|A v_k\|} \]
Le quotient de Rayleigh
Pour estimer la valeur propre à partir du vecteur \(v\), on utilise le quotient de Rayleigh que nous ne développerons pas ici (c’est une notion hors programme). Retiens simplement qu’il s’agit d’un produit scalaire qui nous permet d’estimer notre valeur propre dominante.
Pour information, Si \(\|v\| = 1\), il se simplifie en :
\[ \lambda \approx \langle Av, v \rangle = v^T A v \] (Formule que tu auras peut-être déjà croisée en exercice !)
NB : Pour cet algorithme, nous n’utiliserons que des vecteurs \(v\) de norme 1. Inutile donc de connaître la formule générale du quotient de Rayleigh, retiens simplement la formule dans le cas particulier de norme unitaire.
Implémentation en Python
Voici un script utilisant la bibliothèque numpy afin de déterminer une estimation de la valeur propre dominante et qui prend en entrée une matrice \(A\).
Regarde-le, tu trouveras ci-dessous une explication ligne par ligne de ce qui s’y passe !
Voici une analyse détaillée de ce script
Préparation du cœur de l’algorithme
D’abord, on définit la fonction qui prend trois arguments : la matrice \(A\), une tolérance (pour savoir quand s’arrêter) et un nombre maximum d’itérations (pour éviter que le programme ne boucle à l’infini si la matrice ne converge pas).
On récupère ensuite la dimension de la matrice carrée. Si \(A\) est une matrice \(3 \times 3\), alors \(n = 3\). Cela sert à créer un vecteur de la bonne taille.
On initialise alors un vecteur \(v\) avec des valeurs aléatoires entre 0 et 1. Il est crucial que ce vecteur de départ ne soit pas « orthogonal » au sous-espace propre dominant, sinon l’algorithme ne fonctionnerait pas. En pratique, l’aléatoire garantit presque à 100 % que ce ne sera pas le cas.
Pour poursuivre, on normalise le vecteur pour que sa norme (sa longueur) soit égale à 1.
On crée une variable pour stocker la valeur propre calculée au tour précédent. Au premier tour, elle vaut 0. Elle servira simplement à mesurer la vitesse de l’algorithme.
On définit alors une boucle de maximum 1 000 itérations (puisque nous avons dit plus haut que \(max-iter = 1000\).
Cœur de l’algorithme
Nous voici (enfin) dans le cœur de notre script Python qui applique l’algorithme : on calcule le produit \(A \times v\).
Mathématiquement, à chaque fois qu’on multiplie par \(A\), le vecteur \(w\) s’aligne de plus en plus vers la direction de la valeur propre dominante.
On poursuit alors en calculant le quotient de Rayleigh : puisque \(w = Av\) et que \(v\) est de norme \(1\), le produit scalaire \(v^T Av\) donne une estimation directe de la valeur propre \(\lambda\) que nous estimons (cf. formule de Rayleigh exposée plus haut).
On normalise le nouveau vecteur \(w\) pour obtenir le « prochain » \(v\). Cette étape est vitale, puisque si la valeur propre est supérieure à \(1\), les nombres deviendraient si grands après quelques itérations que l’ordinateur ne saurait les supporter.
Issue du script
On réalise enfin un test de convergence : on regarde si l’écart entre la valeur actuelle et celle du tour précédent est devenu minuscule (inférieure à \(10^{-9}\)).
Si c’est le cas, on considère qu’on a trouvé la solution que l’on demande à l’algorithme de nous retourner.
Si on n’a pas encore convergé, on met à jour la valeur précédente avec la valeur actuelle et on repart pour un tour de boucle.
On ajoute enfin une dernière ligne qui renvoie le résultat même si on a atteint les \(1000\) itérations sans atteindre la précision demandée, dans le cas où l’algorithme n’aurait pas abouti.
Points de vigilance
Analyse de la convergence
La vitesse de convergence est dictée par le « ratio d’amortissement » suivant :
\[ \rho = \left| \frac{\lambda_2}{\lambda_1} \right| \]
Nous admettrons ici ce résultat. Retiens simplement que plus \(\rho\) est proche de 0, plus la convergence est rapide. À l’inverse, si \(\lambda_1\) et \(\lambda_2\) sont très proches, l’algorithme nécessite un grand nombre d’itérations, ce qui peut être problématique pour des matrices de très grande dimension.
Vecteur initial
Si \(v_0\) est orthogonal à \(u_1\) (\(\alpha_1 = 0\)), l’algorithme échoue. En pratique, l’initialisation aléatoire règle ce problème, puisqu’il est quasiment impossible que ces vecteurs soient orthogonaux, mais il est important de connaître cette condition de réalisation de l’algorithme.
Stabilité numérique
La normalisation est l’étape qui permet d’éviter ce qu’on appelle en informatique l’overflow (= dépassement de capacité mémoire). Concrètement, sans cette étape, le script risque de bugger lorsqu’on l’exécute sur un ordinateur.
Conclusion
En définitive, l’algorithme du calcul des valeurs propres dominantes repose sur un cas particulier du quotient de rayleigh (complètement hors programme). Néanmoins, la formule par laquelle nous passons (cas particulier de la formule générale) est souvent abordée dans les sujets. C’est pourquoi il n’est pas improbable qu’un tel sujet sur les valeurs propres dominantes refasse son apparition, en incluant un calcul en Python d’une valeur approchée.
Tu peux retrouver ici le méga-répertoire qui contient toutes les annales de concours et les corrigés. Tu peux également accéder ici à toutes nos autres ressources mathématiques !
