Décomposition SVD

La décomposition en valeurs singulières d'une matrice A de dimension M×N s'écrit A = U W V^T, où U est une matrice orthogonale de dimension M×M et V une matrice orthogonale de dimension N×N. Les éléments diagonaux de la matrice W sont des nombres non négatifs ordonnés de façon décroissante, tandis que tous les éléments hors diagonale sont nuls.

Comme la matrice W est principalement composée de zéros, seules les min(M,N) premières colonnes de la matrice U sont nécessaires pour obtenir la matrice A. De même, seules les min(M,N) premières lignes de la matrice V^T interviennent dans le produit. Ces colonnes et lignes sont appelées respectivement vecteurs singuliers gauche et droit.

La décomposition en valeurs singulières possède de nombreuses propriétés utiles. Elle permet, par exemple, de résoudre des systèmes d'équations linéaires sous-déterminés et surdéterminés, d'effectuer l'inversion et la pseudo-inversion de matrices, de calculer le conditionnement d'une matrice, d'orthogonaliser des systèmes vectoriels et de calculer le complément orthogonal.

Revue des algorithmes existants

L'algorithme de Jacobi fut l'un des premiers à effectuer la décomposition en valeurs singulières (SVD). Il réduit une matrice rectangulaire à une matrice diagonale par une suite de rotations élémentaires. Cette méthode permet de trouver toutes les valeurs singulières avec une grande précision, même les plus petites. Cependant, ses performances étant relativement faibles, les algorithmes QR itératifs ont gagné en popularité.

Les algorithmes SVD modernes les plus répandus reposent sur la réduction de la matrice à une forme bidiagonale par transformation orthogonale (ce problème est suffisamment simple et sa résolution ne requiert qu'un nombre fini d'opérations) et sa diagonalisation par un algorithme QR itératif. Généralement, les algorithmes diffèrent uniquement par leur implémentation de l'algorithme QR. Initialement, la famille d'algorithmes basée sur l'algorithme de Golub-Kahan-Reinsch était la plus répandue. C'est cette même méthode qui a été implémentée dans LINPACK/EISPACK. Ses avantages résident dans sa simplicité et sa concision (300 lignes de code). Malheureusement, elle ne présente pas d'autres atouts. Bien que cet algorithme soit capable de résoudre les problèmes, sa convergence et sa précision sont relativement faibles.

La bibliothèque LAPACK contient deux algorithmes de décomposition en valeurs singulières corrigeant les défauts de l'algorithme LINPACK. Le premier (sous-programmes xBDSQR et xGESVD), prototype de l'algorithme décrit ici, offre une meilleure précision et une convergence plus rapide que son équivalent LINPACK et le remplace donc.

Il convient de noter que ce nouvel algorithme trouve les petites valeurs singulières d'une matrice bidiagonale avec une meilleure précision. La variante LINPACK de l'itération QR trouve toutes les valeurs singulières avec la même erreur absolue. La plus grande valeur singulière est obtenue avec une précision proche de la précision machine. Ceci est suffisant pour les problèmes où l'erreur absolue sur les valeurs singulières est critique : résolution de systèmes d'équations linéaires, inversion de matrices, etc. Cependant, il arrive que de plus petites valeurs singulières soient importantes, même si leur erreur relative s'avère trop importante. Par exemple, si la première valeur singulière est égale à 1 et a été trouvée avec une erreur absolue de 10?¹? (15 chiffres significatifs), la valeur singulière égale à 10^-10, trouvée avec la même erreur absolue, n'aura que 5 chiffres significatifs. Les nouveaux algorithmes trouvent toutes les valeurs avec la même erreur relative, de sorte que les deux valeurs singulières auront 15 chiffres significatifs.

Malheureusement, les erreurs de réduction d'une matrice quelconque à une forme bidiagonale annulent cet avantage : elles corrompent les petites valeurs singulières, qui peuvent être trouvées avec précision grâce au nouvel algorithme. Par conséquent, si votre matrice n'est pas bidiagonale, vous ne bénéficierez du nouvel algorithme que pour son meilleur avantage : une convergence plus rapide. Cependant, certains types de matrices peuvent être réduits à une forme bidiagonale sans modifier les petites valeurs singulières. De tels algorithmes ne sont pas présentés sur le site pour le moment. Si la précision des petites valeurs singulières est essentielle, il est conseillé de rechercher des ressources théoriques et de les implémenter soi-même. Il existe peu d'algorithmes conçus dans ce domaine (espérons que cela s'améliorera).

Comme indiqué précédemment, la bibliothèque LAPACK contient deux algorithmes de décomposition en valeurs singulières (SVD). Le second, dit «diviser pour régner», divise la tâche de décomposition SVD d'une grande matrice bidiagonale en plusieurs sous-tâches résolues par l'algorithme QR. Cet algorithme offre de meilleures performances que l'algorithme QR pour les grandes matrices. Par exemple, la décomposition en valeurs singulières d'une matrice carrée par «diviser pour régner» est 2 à 4 fois plus rapide pour N = 100 qu'avec un simple algorithme QR (temps de réduction de la matrice à une forme bidiagonale inclus), et 6 à 7 fois plus rapide pour N = 1000. Cependant, sa complexité technique rend son intégration à la bibliothèque ALGLIB peu probable dans un avenir proche. Si la performance de la décomposition en valeurs singulières de grandes matrices est essentielle à vos applications, nous vous recommandons d'utiliser la bibliothèque LAPACK.

Algorithme

Comme décrit précédemment, les algorithmes modernes de décomposition singulière réduisent la matrice à une forme bidiagonale, puis la diagonalisent à l'aide de l'algorithme QR. Cette méthode simple est assez efficace, mais il est possible d'accélérer considérablement l'algorithme. Le schéma d'algorithme amélioré décrit ci-dessous est tiré de la bibliothèque LAPACK (sous-programme xGESVD).

Cet algorithme fonctionne bien pour les matrices carrées et quasi-carrées. Si la matrice source est rectangulaire et allongée, au lieu de la réduire à une forme bidiagonale, on peut utiliser la décomposition LQ ou QR (selon la dimension de la matrice la plus grande) pour la représenter comme le produit d'une matrice orthogonale rectangulaire Q et d'une matrice carrée (triangulaire supérieure ou inférieure). On diagonalise ensuite la matrice carrée ainsi obtenue, beaucoup plus petite que la matrice source. Selon le rapport lignes/colonnes, on obtient un algorithme 2 à 6 fois plus rapide que l'algorithme générique.

La seconde amélioration (plus technique) consiste à optimiser le calcul de la matrice U. Lors du calcul de la matrice U, des opérations sur les colonnes sont effectuées. Ces méthodes ne sont pas efficaces (comparativement à l'utilisation du cache du processeur) lorsqu'on travaille avec des méthodes classiques d'entreposage des matrices par lignes, ce qui est courant dans la plupart des langages de programmation. Si de la mémoire supplémentaire est disponible, on peut calculer la matrice U^T en effectuant des opérations sur les lignes de la matrice. Cette méthode est particulièrement efficace si la matrice U est grande. On peut ensuite transposer la matrice U^T.

Ces deux améliorations nécessitent de la mémoire supplémentaire. S'il est impossible d'allouer suffisamment de mémoire, le programmeur peut interdire l'utilisation de la mémoire supplémentaire, mais cela ralentira considérablement l'algorithme.

Description de la sous-routine

• La sous-routine RMatrixSVD effectue la décomposition SVD d'une matrice rectangulaire M×N. Elle renvoie le tableau des valeurs singulières par ordre décroissant et, le cas échéant, les matrices U et V^T. Dans ce cas, elle peut renvoyer soit uniquement les vecteurs singuliers à gauche, soit uniquement ceux à droite, ainsi que des matrices M×M ou N×N complètes (selon les valeurs des paramètres UNeeded et VTNeeded). Notez que VT contient la matrice V^T.
• Les paramètres de la sous-routine sont décrits plus en détail dans les commentaires du programme.
• Cet algorithme est issu de la bibliothèque LAPACK.