LOBPCG

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.

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

${\displaystyle \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 ${\displaystyle 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

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

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

${\displaystyle w:=Tr,}$

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

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

ou, alternativement de manière moins concise,

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

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

${\displaystyle 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 ${\displaystyle \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,

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

(ou ${\displaystyle \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 ${\displaystyle y}$ appartenant à l'espace vectoriel de dimension deux engendré par les deux vecteurs ${\displaystyle x_{i}}$et ${\displaystyle w_{i}}$.

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 (${\displaystyle x^{i+1}}$devient une fonction de ${\displaystyle x^{i}}$, ${\displaystyle w^{i}}$et ${\displaystyle x^{i-1}}$)

${\displaystyle 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 [2] mais augmente les coûts de calcul, et n'est donc généralement pas recommandé.

Au fur et à mesure que les itérations convergent, les vecteurs ${\displaystyle x^{i}}$ et ${\displaystyle 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 ${\displaystyle x^{i-1}}$ avec un vecteur ${\displaystyle p^{i}}$, qui peut être plus éloigné de ${\displaystyle x^{i}}$, dans la base du sous-espace tridimensionnel ${\displaystyle 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[2]. 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.

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. [2] Même dans le cas trivial ${\displaystyle T=I}$ et ${\displaystyle B=I}$ l'approximation résultante avec ${\displaystyle 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 .

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 .

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é ${\displaystyle (v_{1},\lambda _{1})}$ et enchaîner avec une déflation orthogonale pour trouver les couples suivant. Si on appelle ${\displaystyle 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 ${\displaystyle 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 ${\displaystyle A}$ sauf ${\displaystyle (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 [2], 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.

L'approche par blocs dans LOBPCG remplace les vecteurs uniques ${\displaystyle x^{i},\,w^{i},}$ et ${\displaystyle p^{i}}$ avec des vecteurs-blocs, c'est-à-dire des matrices ${\displaystyle X^{i},\,W^{i},}$ et ${\displaystyle P^{i}}$, où, par exemple, chaque colonne de ${\displaystyle 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 ${\displaystyle X^{i+1}}$ est déterminé par la méthode Rayleigh-Ritz sur le sous-espace couvert par toutes les colonnes de matrices ${\displaystyle X^{i},\,W^{i},}$ et ${\displaystyle P^{i}}$ . Chaque colonne de ${\displaystyle W^{i}}$ est calculée simplement comme le résidu préconditionné pour chaque colonne de ${\displaystyle X^{i}.}$ La matrice ${\displaystyle P^{i}}$ est déterminée de telle sorte que les sous-espaces engendrés par les colonnes de ${\displaystyle [X^{i},\,P^{i}]}$ et de ${\displaystyle [X^{i},\,X^{i-1}]}$ sont identiques.

La méthode Rayleigh-Ritz est appliquée sur le sous-espace engendré par toutes les colonnes de matrices ${\displaystyle X^{i},\,W^{i},}$ et ${\displaystyle 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[3] - [4], utilise une factorisation de Cholesky instable mais efficace de la matrice normale, qui est effectuée uniquement sur les matrices individuelles ${\displaystyle W^{i}}$ et ${\displaystyle 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 ${\displaystyle 10^{3}-10^{4}}$, néanmoins le temps de calcul consacré aux orthogonalisations et à la méthode Rayleigh-Ritz commence à dominer.

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[3] - [4] - [2] - [5], 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.

Les progiciels scikit-learn et Megaman[23] utilisent LOBPCG pour adapter le regroupement spectral[24] et l'apprentissage automatique en grande dimension[25] pour les problèmes de réduction de la dimension. NVIDIA a implanté[26] la méthode LOBPCG dans sa bibliothèque nvGRAPH introduite dans CUDA 8. le logiciel Sphynx[27], 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.

LOBPCG est implémenté dans ABINIT[28] (y compris la version CUDA ) et Octopus[29]. 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[30] - [31]. 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.[32] Il existe des versions MATLAB[33] et Julia[34] - [35] - [36] 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[37], Platypus‐QM[38], MFDn[39], ACE-Molecule[40], LACONIC.

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é[41] dans LS-DYNA .

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).

Un filtre passe-bas approché basé sur les itérations de la méthode LOBPCG peut être utilisé pour le débruitage ; voir[42], par exemple, en traitement 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[43] et actuellement testée dans[44] et[45]. Cette dernière approche a ensuite été implémentée dans Python Scikit-learn[46] 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.