LOBPCG

Le Gradient conjugué préconditionné par bloc localement optimal (acronyme anglophone LOBPCG) est une méthode utilisée pour trouver les valeurs propres les plus grandes (ou les plus petites) et les vecteurs propres correspondants d'un problème de valeurs propres généralisé symétrique

$Ax=\lambda Bx,$

pour une paire $(A,B)$ de matrices hermitiennes complexes ou réelles symétriques, où la matrice $B$ est également supposée définie positive.

Caractéristiques principales

• La méthode est dite matrix-free c'est-à-dire qu'elle ne nécessite pas de stocker explicitement les coefficients de la matrice, mais on peut accéder à la matrice en évaluant les produits matrice-vecteur.
• La méthode ne nécessite pas de factorisation de matrice, même pour un problème aux valeurs propres généralisé.
• Le coût de calcul par itération ainsi que l'occupation mémoire sont compétitifs par rapport à la méthode de Lanczos, calculant une seule paire vecteur propre, valeur propre d'une matrice symétrique.
• La convergence linéaire est théoriquement garantie et pratiquement observée.
• On observe une convergence accélérée due au préconditionnement direct, contrairement à la méthode de Lanczos, la matrice du préconditionneur peut être fixe ou varier d'une itération à l'autre mais rester symétrique et définie positive.
• La méthode peut tirer parti de l'utilisation de techniques de décomposition de domaine et d'un préconditionneur multigrille.
• La méthode peut être redémarrée à chaud (c'est-à-dire que l'on utilise des vecteurs propres approchés venant d'une autre méthode pour initier le calcul) et elle calcule une approximation des vecteurs propres, valeurs propres à chaque itération.
• La méthode est plus stable numériquement que la méthode de Lanczos, et peut fonctionner en arithmétique de faible précision.
• Elle est facile à mettre en œuvre, avec de nombreuses implantations existantes.
• La méthode repose sur la notion de bloc et permet ainsi d'utiliser des opérations matrice-matrice très efficaces, par exemple BLAS 3.
• La taille des blocs peut être ajustée pour équilibrer la vitesse de convergence par rapport aux coûts de calcul des orthogonalisations et au coût de la méthode Rayleigh-Ritz à chaque itération.

Algorithme

Version à un vecteur

On commence par décrire un algorithme qui recherché un seul couple vecteur propre / valeur propre, par exemple celui correspondant à la plus petite valeur propre.

Préliminaires : rappel sur la méthode de montée ou descente de gradient pour les problèmes de valeurs propres

La méthode effectue une maximisation (ou minimisation) itérative du quotient de Rayleigh généralisé

$\rho (x):=\rho (A,B;x):={\frac {x^{T}Ax}{x^{T}Bx}},$

ce qui permet de trouver les paires vecteurs propres, valeurs propres les plus grandes (ou les plus petites) de $Ax=\lambda Bx.$

La direction utilisée par l'algorithme du gradient est précisément le gradient du quotient de Rayleigh généralisé, direction proportionnelle au vecteur

$r:=Ax-\rho (x)Bx,$

appelé vecteur propre résiduel. Si un préconditionneur $T$ est disponible, il est appliqué au résidu et donne le vecteur

$w:=Tr,$

appelé résidu préconditionné. Sans préconditionnement, on pose $T:=I$ et donc $w:=r$. La méthode itérative s'écrit

$x^{i+1}:=x^{i}+\alpha ^{i}T(Ax^{i}-\rho (x^{i})Bx^{i}),$

ou, alternativement de manière moins concise,

$x^{i+1}:=x^{i}+\alpha ^{i}w^{i},\,$

$w^{i}:=Tr^{i},\,$

$r^{i}:=Ax^{i}-\rho (x^{i})Bx^{i},$

Cette méthode est connue sous le nom de montée (ou descente) de gradient préconditionnée, où le scalaire $\alpha ^{i}$ est le pas d'itération. La taille optimale du pas peut être déterminée en maximisant le quotient de Rayleigh, c'est-à-dire,

$x^{i+1}:=\arg \max _{y\in span\{x^{i},w^{i}\}}\rho (y)$

(ou $\arg \min$ en cas de minimisation), auquel cas la méthode est dite localement optimale.

On note que le maximum (ou minimum) est pris pour $y$ appartenant à l'espace vectoriel de dimension deux engendré par les deux vecteurs $x_{i}$et $w_{i}$.

Formule de récurrence à trois termes

Pour accélérer considérablement la convergence de la montée (ou descente) de gradient préconditionnée localement optimale, un vecteur supplémentaire peut être ajouté à la formule de récurrence à deux termes pour la rendre à trois termes ($x^{i+1}$devient une fonction de $x^{i}$, $w^{i}$et $x^{i-1}$)

$x^{i+1}:=\arg \max _{y\in span\{x^{i},w^{i},x^{i-1}\}}\rho (y)$

La maximisation/minimisation du quotient de Rayleigh se fait maintenant dans un espace vectoriel de dimension trois et peut être effectuée numériquement par la méthode Rayleigh-Ritz. L'ajout de vecteurs supplémentaires, voir, par exemple, l'extrapolation de Richardson, n'entraîne pas d'accélération significative mais augmente les coûts de calcul, et n'est donc généralement pas recommandé.

Améliorations de la stabilité numérique

Au fur et à mesure que les itérations convergent, les vecteurs $x^{i}$ et $x^{i-1}$ deviennent presque linéairement dépendants, ce qui entraîne une perte de précision et rend la méthode Rayleigh-Ritz numériquement instable en présence d'erreurs d'arrondi. La perte de précision peut être évitée en substituant le vecteur $x^{i-1}$ avec un vecteur $p^{i}$, qui peut être plus éloigné de $x^{i}$, dans la base du sous-espace tridimensionnel $span\{x^{i},w^{i},x^{i-1}\}$, tout en gardant le sous-espace inchangé et en évitant l'orthonormalisation ou toute autre opération supplémentaire. De plus, l'orthonormalisation de la base du sous-espace tridimensionnel peut être nécessaire pour les problèmes de valeurs propres mal conditionnés afin d'améliorer la stabilité et la précision atteignable.

Analogues du sous-espace de Krylov

Il s'agit d'une version à vecteur unique de la méthode LOBPCG - une possible généralisation des solveurs linéaires à gradient conjugué préconditionnés au cas des problèmes de valeurs propres symétriques. Même dans le cas trivial $T=I$ et $B=I$ l'approximation résultante avec $i>3$ sera différent de celui obtenu par l'algorithme de Lanczos, bien que les deux approximations appartiendront au même sous-espace de Krylov.

Scénarios d'utilisation pratiques

L'extrême simplicité et la haute efficacité de la version à vecteur unique de LOBPCG la rendent attrayante pour les applications de calcul de valeurs propres en environnement contraint en ressources de calcul, par exemple pour la détection d'anomalies en temps réel basée sur le partitionnement spectral. On utilise aussi pour le partitionnement de graphe sur ASIC ou FPGA, ou encore pour la modélisation de phénomènes physiques sur les supercalculateurs exascale du classement TOP500.

version par bloc de l'algorithme

Résumé

On cherche maintenant à calculer plusieurs couples propres simultanément. On pourrait appliquer la méthode LOBPCG à un seul vecteur (décrite ci-dessus) pour avoir le couple correspondant à la plus grande valeur propre (en module), noté $(v_{1},\lambda _{1})$ et enchaîner avec une déflation orthogonale pour trouver les couples suivant. Si on appelle $A$ la matrice (supposée symétrique, mais on peut généraliser on cas non-symétrique) dont on cherche les couples propres, rappelons que la déflation orthogonale consiste à introduire une matrice $A'=A-{\frac {\lambda _{1}}{|v_{1}|^{2}}}v_{1}^{T}v_{1}$dont on peut montrer qu'elle a les mêmes couples propres que $A$ sauf $(v_{1},\lambda _{1})$. Il suffit d'appliquer à nouveau la méthode LOBPCG à un seul vecteur pour trouver un deuxième couple propre et de répéter cette opération autant de fois que nécessaire.

Cependant, avec approche, chaque nouveau couple propre calculé de manière approchée dépend du calcul de tous les autres couples propres précédents et ainsi les erreurs de calcul s'accumulent et les estimations des couples propres deviennent de moins en moins bonnes. Pour améliorer cette situation, la méthode LOBPCG par bloc se propose de calculer plusieurs couples propres ensemble d'un seul bloc, de manière rapide, précise et robuste des vecteurs propres. Elle améliore également les cas où l'on a plusieurs plusieurs couples propres très proches les uns des autres qui seraient très mal estimés par la déflation orthogonale (convergence très lente). La taille du bloc peut être ajustée pour équilibrer la stabilité numérique par rapport à la vitesse de convergence et aux coûts de calcul des orthonormalisations et de la méthode Rayleigh-Ritz à chaque itération.

Principe de la méthode

L'approche par blocs dans LOBPCG remplace les vecteurs uniques $x^{i},\,w^{i},$ et $p^{i}$ avec des vecteurs-blocs, c'est-à-dire des matrices $X^{i},\,W^{i},$ et $P^{i}$, où, par exemple, chaque colonne de $X^{i}$ se rapproche de l'un des vecteurs propres. Toutes les colonnes sont itérées simultanément, et la prochaine matrice de vecteurs propres approximatifs $X^{i+1}$ est déterminé par la méthode Rayleigh-Ritz sur le sous-espace couvert par toutes les colonnes de matrices $X^{i},\,W^{i},$ et $P^{i}$. Chaque colonne de $W^{i}$ est calculée simplement comme le résidu préconditionné pour chaque colonne de $X^{i}.$ La matrice $P^{i}$ est déterminée de telle sorte que les sous-espaces engendrés par les colonnes de $[X^{i},\,P^{i}]$ et de $[X^{i},\,X^{i-1}]$ sont identiques.

Compromis entre stabilité numérique et efficacité

La méthode Rayleigh-Ritz est appliquée sur le sous-espace engendré par toutes les colonnes de matrices $X^{i},\,W^{i},$ et $P^{i}$, où en principe une base peut être choisie arbitrairement. Cependant, la méthode Rayleigh-Ritz peut en pratique devenir numériquement instable si certains des vecteurs de base ne sont qu'approximativement linéairement indépendant.

Le choix de la base de sous-espace devient crucial autant pour d'assurer la stabilité numérique de la méthode Rayleigh-Ritz que pour réduire à des coûts de calcul. L'approche sans doute la plus stable consiste à rendre les vecteurs de base orthogonaux, par exemple, par le processus d'orthonormalisation de Gram-Schmidt, opération coûteuse en temps de calcul. Par exemple, les implémentations LOBPCG, utilise une factorisation de Cholesky instable mais efficace de la matrice normale, qui est effectuée uniquement sur les matrices individuelles $W^{i}$ et $P^{i}$, plutôt que sur l'ensemble du sous-espace. Sur une architecture de calcul moderne, la quantité de mémoire disponible est suffisante pour se permettre d'utiliser des tailles de bloc de l'ordre de $10^{3}-10^{4}$, néanmoins le temps de calcul consacré aux orthogonalisations et à la méthode Rayleigh-Ritz commence à dominer.

Verrouillage des vecteurs propres précédemment convergés

Quand on utilise une méthode par bloc itérative pour approcher les vecteurs propres, il arrive très souvent qu'à l'intérieur d'un bloc, certains vecteurs propres convergent plus rapidement que d'autres. Le verrouillage des vecteurs propres déjà convergés consiste à les retirer de la boucle itérative, afin d'éliminer des calculs inutiles et d'améliorer la stabilité numérique. Une simple suppression d'un vecteur propre peut probablement entraîner la formation de son double dans des vecteurs toujours itératifs. Le fait que les vecteurs propres d'une matrice symétrique soient orthogonaux par paire suggère de garder tous les vecteurs itératifs orthogonaux aux vecteurs verrouillés.

Le verrouillage peut être mis en œuvre différemment en maintenant la précision et la stabilité numériques tout en minimisant les coûts de calcul. Par exemple, les implémentations LOBPCG, séparent le verrouillage dur, c'est-à-dire une déflation par restriction, où les vecteurs propres verrouillés servent d'entrée de code et ne changent pas, du verrouillage doux, où les vecteurs verrouillés ne participent pas à l'étape itérative généralement la plus coûteuse en calcul des résidus, mais cependant, participent pleinement à la méthode Rayleigh-Ritz et peuvent donc être modifiés par la méthode Rayleigh-Ritz.

Théorie et pratique de la convergence

La méthode LOBPCG, par construction, assure une minimisation du quotient de Rayleigh au moins aussi rapide que la descente de gradient appliqué à un du bloc, dont on connaît la convergence théorique. Chaque vecteur propre est un point stationnaire du quotient de Rayleigh, où le gradient s'annule. Ainsi, la descente du gradient peut ralentir à proximité de n'importe quel vecteur propre, cependant, soit la méthode converge vers le vecteur propre avec un taux de convergence linéaire ou, si ce vecteur propre est un point col, le quotient de Rayleigh itératif est susceptible de converger vers la deuxième plus basse valeur propre. La convergence dans le pire cas a été déterminée et dépend de l'écart relatif entre la valeur propre et le reste du spectre de la matrice et de la qualité du préconditionneur, s'il est présent.

Pour une matrice quelconque, il n'y a évidemment aucun moyen de prédire les vecteurs propres et donc de générer les premières approximations qui fonctionnent toujours bien. La solution itérative de LOBPCG peut être sensible aux valeurs initiales des vecteurs propres, et par exemple, mettre plus de temps à converger et ralentir lorsque les paires propres suivantes doivent être calculées. De plus, en théorie, on ne peut pas nécessairement garantir la convergence vers le plus petit couple propre, bien que la probabilité de l'échec soit nulle. Une fonction gaussienne aléatoire de bonne qualité avec une moyenne nulle est généralement prise pour générer les approximations initiales de la méthode LOBPCG. Pour fixer les valeurs approchées initiales, on peut choisir une graine fixée pour le générateur de nombres aléatoires.

Contrairement à la méthode de Lanczos, LOBPCG présente rarement une convergence superlinéaire asymptotique dans la pratique.

Application à l'analyse en composantes principales (ACP) et à la décomposition en valeur singulière (SVD)

la méthode LOBPCG peut être facilement adaptée au calcul de plusieurs plus grandes valeurs singulières et les vecteurs singuliers correspondants et ainsi être utilisée pour l'analyse en composante principale (ACP). Elle est implantée dans le module Python SciPy depuis la version 1.4.0.

Implémentations logicielles

Andrew Knyazev, l'auteur original de la méthode LOBPCG, a publié une implantation logicielle de référence appelée Block Locally Optimal Eigenvalue Xolvers (BLOPEX) avec des interfaces vers les bibliothèques d'algèbre linéaire parallèle PETSc, hypre et Parallel Hierarchical Adaptive MultiLevel method (PHAML). D'autres implantations sont disponibles dans les logiciels GNU Octave, MATLAB, Java, Anasazi (Trilinos), SLEPc, SciPy, Julia, MAGMA, PyTorch, Rust, OpenMP et OpenACC, CuPy (Une bibliothèque de tableaux compatible NumPy pour les processeurs graphiques), et NVIDIA AMGX. Il existe également une implantation de LOBPCG compatible avec TensorFlow,

Applications

Exploration de données

Les progiciels scikit-learn et Megaman utilisent LOBPCG pour adapter le regroupement spectral et l'apprentissage automatique en grande dimension pour les problèmes de réduction de la dimension. NVIDIA a implanté la méthode LOBPCG dans sa bibliothèque nvGRAPH introduite dans CUDA 8. le logiciel Sphynx, implante un algorithme de partitionnement de graphe parallèle hybride (CPU et GPU) à mémoire distribuée et partagée, il utilise le regroupement spectral pour le partitionnement de graphe, en calculant les vecteurs propres sur la matrice laplacienne du graphe en utilisant la méthode LOBPCG du package Anasazi de la bibliothèques Trilinos.

Sciences des matériaux

LOBPCG est implémenté dans ABINIT (y compris la version CUDA) et Octopus. Il a été utilisé pour des matrices de plusieurs milliards de dollars par les finalistes du prix Gordon Bell, sur le supercalculateur Earth Simulator au Japon. Le modèle Hubbard pour les systèmes électroniques fortement corrélés pour comprendre le mécanisme derrière la supraconductivité utilise LOBPCG pour calculer l'état fondamental de l'hamiltonien sur l'ordinateur K. Il existe des versions MATLAB et Julia de LOBPCG pour les équations de Kohn-Sham et la théorie de la fonctionnelle de la densité (DFT) utilisant la base des ondes simples. Les implémentations récentes incluent TTPY, Platypus‐QM, MFDn, ACE-Molecule, LACONIC.

Mécaniques des fluides

L'implantation de LOBPCG de la bibliothèque BLOPEX est utilisée pour la configuration du préconditionneur dans la bibliothèque de solveurs BDDCML (Multilevel Balancing Domain Decomposition by Constraints), qui est incorporé dans OpenFTL (Open Finite Element Template Library) et utilisé par le simulateur Flow123d pour étudier les écoulements d'eau souterraine, de solutés et de transport de chaleur dans les milieux poreux fracturés. LOBPCG a été implémenté dans LS-DYNA.

Electromagnétisme

LOBPCG est l'un des principales méthodes de recherche de valeurs propres dans PYFEMax et dans le logiciel d'éléments finis multiphysique haute performance Netgen/NGSolve. L'implantation de LOBPCG fournie par la bibliothèque hypre est intégré à la bibliothèque C++ MFEM (simulation par la méthode des éléments finis), qui est utilisée dans de nombreux projets, notamment BLAST, XBraid, VisIt, xSDK, et par l'institut FASTMath du SciDAC pour la conception de discrétisations exascales efficaces (CEED).

Débruitage

Un filtre passe-bas approché basé sur les itérations de la méthode LOBPCG peut être utilisé pour le débruitage ; voir, par exemple, en traitement d'images.

Segmentation d'images

La segmentation d'images par regroupement spectral effectue un plongement de faible dimension à l'aide d'une matrice d'affinité entre les pixels, suivi d'un regroupement des composants des vecteurs propres dans l'espace de faible dimension, par exemple en utilisant le graphe Laplacien pour le filtrage bilatéral. La segmentation d'images via un partitionnement de graphe spectral par LOBPCG avec préconditionnement multigrille a été proposée pour la première fois dans et actuellement testée dans. Cette dernière approche a ensuite été implémentée dans Python Scikit-learn qui utilise l'implantation de LOBPCG fournie par SciPy avec un préconditionnement algébrique multigrille pour résoudre le problème des valeurs propres pour le graphe Laplacien.

Liens externes

• LOBPCG dans MATLAB
• LOBPCG à Octave
• LOBPCG dans SciPy
• LOBPCG en Java chez Google Code
• LOBPCG dans Block Locally Optimal Eigenvalue Xolvers (BLOPEX) sur GitHub et archivé sur Google Code