Algèbre linéaire numérique

L'algèbre linéaire numérique a été développée par des mathématiciens pionniers de l'informatique comme John von Neumann, Alan Turing, James H. Wilkinson, Alston Scott Householder, George Forsythe et Heinz Rutishauser, afin d'utiliser les premiers ordinateurs dans la résolution de problèmes de mathématiques continues, comme les problèmes de balistique et les systèmes d'équations aux dérivées partielles^[2]. La première tentative sérieuse de minimiser les erreurs informatiques lors de l'application d'algorithmes à des données réelles est celle de John von Neumann et Herman Goldstine en 1947^[3]. Ce domaine s'est développé à mesure que la technologie a permis aux chercheurs de résoudre des problèmes complexes sur des matrices de très grande taille et de haute précision, et certains algorithmes numériques ont acquis une importance croissante grâce à des technologies comme le calcul parallèle qui en ont fait des approches pratiques pour la résolution de problèmes scientifiques^[2].

Pour de nombreux problèmes d'algèbre linéaire appliquée, il est utile de considérer une matrice comme une concaténation de vecteurs colonnes. Par exemple, lors de la résolution du système linéaire ${\displaystyle x=A^{-1}b}$, plutôt que de comprendre x comme le produit de la matrice inverse ${\displaystyle A^{-1}}$ avec le vecteur b, il est utile de considérer x comme le vecteur des coefficients du développement linéaire de b dans la base formée par les colonnes de A^[1]:8. Considérer les matrices comme une concaténation de colonnes est également une approche pratique pour les algorithmes matriciels. En effet, ces algorithmes contiennent fréquemment deux boucles imbriquées : l’une parcourant les colonnes d’une matrice A, et l’autre ses lignes. Par exemple, pour les matrices ${\displaystyle A^{m\times n}}$ et les vecteurs ${\displaystyle x^{n\times 1}}$ et ${\displaystyle y^{m\times 1}}$, on peut utiliser la perspective de partitionnement en colonnes pour calculer y := Ax + y comme

for q = 1:n
 for p = 1:m
  y(p) = A(p,q)*x(q) + y(p)
 end
end

La décomposition en valeurs singulières d'une matrice ${\displaystyle A^{m\times n}}$ consiste à trouver la forme ${\displaystyle A=U\Sigma V^{\ast }}$ où U et V sont unitaires, et ${\displaystyle \Sigma }$ est diagonale. Les éléments diagonaux de ${\displaystyle \Sigma }$ sont appelées les valeurs singulières de A. Parce que les valeurs singulières sont les racines carrées des valeurs propres de ${\displaystyle AA^{\ast }}$, il existe un lien étroit entre la décomposition en valeurs singulières et la décomposition en valeurs propres. Cela signifie que la plupart des méthodes de calcul de la décomposition en valeurs singulières sont similaires aux méthodes de calcul des valeurs propres^[1]:36. La méthode peut-être la plus courante utilise les transformations de Householder^[1]:253.

La factorisation QR d'une matrice ${\displaystyle A^{m\times n}}$ est le produit d'une matrice ${\displaystyle Q^{m\times m}}$ et d'une matrice ${\displaystyle R^{m\times n}}$ de sorte que A = QR, où Q est orthogonal et R est triangulaire supérieure^{[1]:50,[4]:223}. Les deux principaux algorithmes de calcul des factorisations QR sont le procédé de Gram-Schmidt et la transformation de Householder. La factorisation QR est souvent utilisée pour résoudre les problèmes de moindres carrés linéaires et les problèmes de valeurs propres (via l'algorithme QR itératif).

La factorisation LU d'une matrice A consiste à construire une matrice triangulaire inférieure L et une matrice triangulaire supérieure U, telles que A = LU. La matrice U est obtenue par une procédure de triangularisation supérieure qui consiste à multiplier à gauche A par une série de matrices ${\displaystyle M_{1},\ldots ,M_{n-1}}$ pour former le produit ${\displaystyle M_{n-1}\cdots M_{1}A=U}$, de sorte que de manière équivalente ${\displaystyle L=M_{1}^{-1}\cdots M_{n-1}^{-1}}$^{[1]:147,[4]:96}.

La diagonalisation d'une matrice ${\displaystyle A^{m\times m}}$ est ${\displaystyle A=X\Lambda X^{-1}}$, où les colonnes de X sont les vecteurs propres de A, et ${\displaystyle \Lambda }$ est une matrice diagonale dont les éléments diagonaux sont les valeurs propres correspondantes de A^[1]:33. Il n'existe pas de méthode directe pour trouver la décomposition en valeurs propres d'une matrice quelconque. Comme il est impossible d'écrire un programme qui trouve les racines exactes d'un polynôme quelconque en temps fini, tout solveur de valeurs propres général doit nécessairement être itératif^[1]:192.

Du point de vue de l'algèbre linéaire numérique, l'élimination de Gauss est une procédure permettant de factoriser une matrice A en sa factorisation LU, ce que l'élimination de Gauss accomplit en multipliant à gauche A par une suite de matrices ${\displaystyle L_{m-1}\cdots L_{2}L_{1}A=U}$ jusqu'à ce que U soit triangulaire supérieure et L triangulaire inférieure, où ${\displaystyle L\equiv L_{1}^{-1}L_{2}^{-1}\cdots L_{m-1}^{-1}}$^{[1] :148}. Les programmes naïfs d'élimination gaussienne sont connus pour leur grande instabilité et produisent des erreurs considérables lorsqu'ils sont appliqués à des matrices comportant de nombreux chiffres significatifs^[2]. La solution la plus simple consiste à introduire le pivot, ce qui permet d'obtenir un algorithme d'élimination de Gauss modifié et stable^{[1] :151}.

L'algèbre linéaire numérique considère généralement les matrices comme une concaténation de vecteurs colonnes. Afin de résoudre le système linéaire ${\displaystyle x=A^{-1}b}$, l'approche algébrique traditionnelle consiste à considérer x comme le produit de ${\displaystyle A^{-1}}$ avec b. L'algèbre linéaire numérique interprète plutôt x comme le vecteur des coefficients de la décomposition de b dans la base formée par les colonnes de A^[1]:8.

De nombreuses décompositions différentes peuvent être utilisées pour résoudre le problème linéaire, en fonction des caractéristiques de la matrice A et des vecteurs x et b, ce qui peut rendre une factorisation beaucoup plus facile à obtenir que d'autres. Si A = QR est une factorisation QR de A, alors de manière équivalente on a ${\displaystyle Rx=Q^{\ast }b}$. Cela se calcule aussi facilement qu'une factorisation matricielle^{[1] :54}. Si ${\displaystyle A=X\Lambda X^{-1}}$ est une décomposition spectrale A, et on cherche à trouver b tel que b = Ax, avec ${\displaystyle b'=X^{-1}b}$ et ${\displaystyle x'=X^{-1}x}$, alors on a ${\displaystyle b'=\Lambda x'}$^{[1] :33}. Ceci est étroitement lié à la résolution du système linéaire par décomposition en valeurs singulières, car les valeurs singulières d'une matrice sont les valeurs absolues de ses valeurs propres, qui sont également équivalentes aux racines carrées des valeurs absolues des valeurs propres de la matrice de Gram ${\displaystyle X^{*}X}$. Et si A = LU est une factorisation LU de A, alors Ax = b peut être résolu par un algorithme de descente-remontée à l'aide des matrices triangulaires Ly = b et Ux = y^{[1] :147,[4]:99}.

Les décompositions matricielles suggèrent plusieurs façons de résoudre le système linéaire r = b − Ax, où l'on cherche à minimiser r, comme dans le problème de régression. L'algorithme QR résout ce problème en calculant la factorisation QR réduite de A et en la réarrangeant pour obtenir ${\displaystyle {\widehat {R}}x={\widehat {Q}}^{\ast }b}$. Ce système triangulaire supérieur peut alors être résolu pour x. La décomposition en valeurs singulières (SVD) propose également un algorithme pour obtenir la solution du problème exprimé comme un problème de moindres carrés linéaires. En calculant la décomposition SVD réduite ${\displaystyle A={\widehat {U}}{\widehat {\Sigma }}V^{\ast }}$ puis en calculant le vecteur ${\displaystyle {\widehat {U}}^{\ast }b}$, on ramène le problème des moindres carrés à un système diagonal simple^[1]:84. Le fait que les factorisations QR et la SVD permettent d'obtenir des solutions par moindres carrés signifie qu'en plus de la méthode classique des équations normales pour résoudre les problèmes de moindres carrés, ces problèmes peuvent également être résolus par des méthodes incluant l'algorithme de Gram-Schmidt et les méthodes de Householder.

On considère un problème sous forme d'une fonction ${\displaystyle f:X\to Y}$, où X est un espace vectoriel normé de données et Y est un espace vectoriel normé de solutions. Pour un certain point de données ${\displaystyle x\in X}$. On dit que le problème est mal conditionné si une petite perturbation de x entraîne une grande variation de la valeur de f ( x ). On peut quantifier cela en définissant un nombre de conditionnement qui mesure la qualité du conditionnement du problème, défini par $${\displaystyle {\widehat {\kappa }}=\lim _{\delta \to 0}\sup _{\|\delta x\|\leq \delta }{\frac {\|\delta f\|}{\|\delta x\|}}.}$$

L'instabilité est la tendance des algorithmes informatiques, dont les calculs reposent sur l'arithmétique à virgule flottante, à produire des résultats très différents de la solution mathématique exacte d'un problème. Lorsqu'une matrice contient des données réelles avec de nombreux chiffres significatifs, de nombreux algorithmes de résolution de problèmes tels que les systèmes linéaires ou l'optimisation par moindres carrés peuvent produire des résultats très imprécis. La création d'algorithmes stables pour les problèmes mal conditionnés est ainsi un enjeu central de l'algèbre linéaire numérique. Par exemple, la stabilité de la triangularisation par transformations de Householder en fait une méthode de résolution particulièrement robuste pour les systèmes linéaires, tandis que l'instabilité de la méthode de la pseudo-inverse pour la résolution des problèmes de moindres carrés justifie le recours à des méthodes de décomposition matricielle, comme la décomposition en valeurs singulières. Certaines méthodes de décomposition matricielle peuvent être instables, mais des modifications simples permettent de les stabiliser ; c'est le cas de la méthode de Gram-Schmidt instable, qui peut être facilement modifiée pour donner la méthode de Gram-Schmidt modifiée stable^[5]:140. Un autre problème classique en algèbre linéaire numérique est la constatation que l'élimination de Gauss est instable, mais devient stable avec l'utilisation d'un pivot partiel.

Les algorithmes itératifs constituent une partie importante de l'algèbre linéaire numérique pour deux raisons principales. Premièrement, de nombreux problèmes numériques importants n'admettent pas de solution directe ; pour trouver les valeurs propres et les vecteurs propres d'une matrice quelconque, seule une approche itérative est possible. Deuxièmement, il existe des algorithmes non itératifs pour une matrice quelconque ${\displaystyle m\times m}$ nécessite un temps en ${\displaystyle O(m^{3})}$, ce qui représente un seuil particulièrement élevé, d'autant que les matrices ne contiennent au plus que ${\displaystyle m^{2}}$ coefficients. Les approches itératives peuvent tirer parti de certaines caractéristiques des matrices pour réduire ce temps de calcul. Par exemple, lorsqu'une matrice est creuse, un algorithme itératif peut s'affranchir de nombreuses étapes qu'une approche directe suivrait nécessairement, même si ces étapes sont redondantes compte tenu de la structure de la matrice.

Le principe de nombreuses méthodes itératives en algèbre linéaire numérique repose sur la projection d'une matrice sur un sous-espace de Krylov de dimension inférieure. Cette projection permet d'approcher les propriétés d'une matrice de grande dimension en calculant itérativement les propriétés équivalentes de matrices similaires, en partant d'un espace de faible dimension et en passant successivement à des dimensions supérieures. Lorsque A est symétrique et que l'on souhaite résoudre le problème linéaire Ax = b, l'approche itérative classique est la méthode du gradient conjugué. Si A n'est pas symétrique, la méthode du résidu minimal généralisé et la méthode du gradient biconjugué sont des exemples de solutions itératives au problème linéaire. Si A est symétrique, l'algorithme de Lanczos permet de résoudre le problème des valeurs propres et des vecteurs propres ; dans le cas général, on utilise plutôt l'algorithme d'Arnoldi.

Plusieurs langages de programmation utilisent des techniques d'optimisation en algèbre linéaire numérique et sont conçus pour implémenter des algorithmes d'algèbre linéaire numérique. Parmi ces langages figurent MATLAB, Analytica, Maple et Mathematica. D'autres langages, non spécifiquement dédiés à l'algèbre linéaire numérique, disposent de bibliothèques fournissant des routines et des fonctions d'optimisation : C et Fortran proposent des packages tels que BASIC (Basic Linear Algebra Subprograms) et LAPACK, Python offre la bibliothèque NumPy et Perl le Perl Data Language. De nombreuses commandes d'algèbre linéaire numérique en R s'appuient sur ces bibliothèques fondamentales, comme LAPACK^[6].