Utilisateur:Mattelbe/Algorithme de multiplication de matrices

De par le fait que la multiplication matricielle est une opération centrale dans de nombreux algorithmes numériques, un travail considérable a été effectué pour rendre les algorithmes de multiplication matricielle efficaces. Les applications de la multiplication matricielle dans des problèmes calculatoires se retrouvent dans de nombreux domaines, notamment le calcul scientifique et la reconnaissance des formes, ainsi que dans des problèmes apparemment sans rapport, tels que le comptage des chemins à travers un graphe. De nombreux algorithmes différents ont été conçus pour multiplier les matrices sur différents types de matériel informatique, y compris les systèmes parallèles et distribués, où les calculs sont répartis sur plusieurs processeurs pouvant être eux-même répartis sur un réseau.

L'application directe de la définition mathématique de la multiplication matricielle donne un algorithme qui prend un temps de l'ordre de $n^{3}$ pour multiplier deux matrices $n\times n$ ($\Theta (n^{3})$ en notation grand O). Bien que de meilleures limites asymptotiques sur le temps nécessaire pour multiplier les matrices sont connues depuis les travaux de Strassen dans les années 1960, on ignore encore la complexité du problème et donc le temps optimal théorique pour ces calculs.

Le résultat de la multiplication matricielle pour deux matrices $A$ et $B$ de tailles respectives $n\times m$ et $m\times p$ est la matrice $C=AB$ de taille $n\times p$ ayant pour éléments

$c_{ij}=\sum _{k=1}^{m}a_{ik}b_{kj}$ .

On peut construire à partir de cette définition un algorithme simple qui boucle sur les indices $i$ de 1 à $n$ et $j$ de 1 à $p$, calculant le résultat précédent à l'aide d'une boucle imbriquée :

• Entrée: matricesA et B
• Soit C une nouvelle matrice de taille appropriée
• Pour i de 1 à n:
  • Pour j de 1 à p:
    • Soit somme = 0
    • Pour k de 1 à m:
      • Affecter somme ← somme + A_ik × B_kj
    • Affecter C_ij ← somme
• Return C

Cet algorithme a une complexité temporelle en $\Theta (nmp)$ (en notation asymptotique).

Une simplification courante pour faciliter l'analyse des algorithmes consiste à supposer que les matrices en entrée sont toutes des matrices carrées de taille $n\times n$. Dans ce cas la complexité temporelle de l'algorithme naïf est en $\Theta \left(n^{3}\right)$, c'est-à-dire cubique.

Comportement du cache

Les trois boucles de la multiplication matricielle itérative peuvent être arbitrairement interchangées sans impacter la validité du résultat ou le temps d'exécution asymptotique. Toutefois, l'ordre peut avoir un impact considérable sur les performances en pratique en raison des modèles d'accès à la mémoire du matériel utilisé et de l'utilisation du cache de l'algorithme ; l'ordre optimal dépend aussi du fait que les matrices sont stockées dans l'ordre des lignes principales, des colonnes principales ou un mélange des deux.

En particulier, dans le cas idéal d'un cache entièrement associatif composé de $M$ octets et $b$ octets par ligne de cache (c'est-à-dire $M/b$ lignes de cache), l'algorithme ci-dessus n'est pas optimal pour $A$ et $B$ stockées dans l'ordre des lignes principales. Lorsque $n>M/b$, chaque itération de la boucle interne (consistant en un balayage simultané à travers une ligne de $A$ et une colonne de $B$) entraîne un échec de cache lors de l'accès à un élément de $B$. Cela signifie que l'algorithme génère $\Theta \left(n^{3}\right)$ échecs de cache dans le pire cas.

Depuis 2010, la vitesse des mémoires par rapport à celle des processeurs est telle que les échecs de cache sont le principal facteur augmentant le temps d'exécution des algorithmes de multiplication de matrices de tailles importantes.

La variante optimale de l'algorithme itératif pour $A$ et $B$ dans la disposition en ligne principale est une version en mosaïque, où la matrice est implicitement divisée en carreaux carrés de taille ${\sqrt {M}}\times {\sqrt {M}}$ :

• Input: matrices A and B
• Let C be a new matrix of the appropriate size
• Pick a tile size T = Θ(Modèle:Radic)
• For I from 1 to n in steps of T:
  • For J from 1 to p in steps of T:
    • For K from 1 to m in steps of T:
      • Multiply A_{I:I+T, K:K+T} and B_{K:K+T, J:J+T} into C_{I:I+T, J:J+T}, that is:
      • For i from I to min(I + T, n):
        • For j from J to min(J + T, p):
          • Let sum = 0
          • For k from K to min(K + T, m):
            • Set sum ← sum + A_ik × B_kj
          • Set C_ij ← C_ij + sum
• Return C

Dans un modèle de cache idéalisé, cet algorithme n'engendre que $\Theta \left({\frac {n^{3}}{b{\sqrt {M}}}}\right)$ échecs de cache.

Le diviseur $b{\sqrt {M}}$ équivaut à plusieurs ordres de grandeur sur les machines modernes, de sorte que les calculs réels dominent le temps de fonctionnement, plutôt que les échecs de cache.

Algorithme diviser pour régner

Une alternative à l'algorithme itératif est un algorithme de type diviser pour régner. Il repose sur le partitionnement des blocs

$C={\begin{pmatrix}C_{11}&C_{12}\\C_{21}&C_{22}\\\end{pmatrix}},\,A={\begin{pmatrix}A_{11}&A_{12}\\A_{21}&A_{22}\\\end{pmatrix}},\,B={\begin{pmatrix}B_{11}&B_{12}\\B_{21}&B_{22}\\\end{pmatrix}}$ ,

qui fonctionne pour toutes les matrices carrées dont les dimensions sont des puissances de deux, c'est-à-dire de forme $2^{n}\times 2^{n}$ pour $n$ donné. Le produit matriciel devient

${\begin{pmatrix}C_{11}&C_{12}\\C_{21}&C_{22}\\\end{pmatrix}}={\begin{pmatrix}A_{11}&A_{12}\\A_{21}&A_{22}\\\end{pmatrix}}{\begin{pmatrix}B_{11}&B_{12}\\B_{21}&B_{22}\\\end{pmatrix}}={\begin{pmatrix}A_{11}B_{11}+A_{12}B_{21}&A_{11}B_{12}+A_{12}B_{22}\\A_{21}B_{11}+A_{22}B_{21}&A_{21}B_{12}+A_{22}B_{22}\\\end{pmatrix}}$

et consiste en huit multiplications de paires de sous-matrices, suivies d'une étape d'addition. L'algorithme diviser pour régner calcule récursivement les plus petites multiplications, en utilisant la multiplication scalaire $c_{11}=a_{11}b_{11}$ comme cas de base.

La complexité de cet algorithme en fonction de $n$ est donnée par la relation de récurrence

$T(1)=\Theta (1)$ ;

$T(n)=8T(n/2)+\Theta (n^{2})$ ,

qui tient compte des huit appels récursifs sur des matrices de taille $n/2$ et $\Theta \left(n^{2}\right)$ pour additionner les quatre paires de matrices résultantes par élément. L'application du théorème sur les récurrences de partition montre que cette récurrence a une solution asymptotique en $\Theta \left(n^{3}\right)$ identique à la complexité de l'algorithme itératif.

Matrices non carrées

Une variante de cet algorithme qui fonctionne pour les matrices de formes arbitraires et qui est plus rapide en pratique divise les matrices en deux au lieu de quatre sous-matrices, comme suit. Fractionner une matrice signifie désormais la diviser en deux parties de taille à peu près égales.

Algorithmes sous-cubiques

Évolution de la puissance ω dans la complexité temporelle en O(n^ω) pour différents algorithmes, en fonction du temps.

Il existe des algorithmes offrant de meilleurs temps d'exécution que les approches naïves. Le premier proposé fut l'algorithme de Strassen, conçu par Volker Strassen en 1969 et souvent appelé « multiplication matricielle rapide ». Il repose sur une manière de multiplier deux matrices 2 × 2 qui ne nécessite que 7 multiplications (au lieu des 8 habituelles), au détriment de plusieurs opérations supplémentaires d'addition et de soustraction. L'application récursive donne un algorithme avec un coût en multiplications de $O(n^{\log _{2}7})\approx O(n^{2.807})$. L'algorithme de Strassen est plus complexe, et sa stabilité numérique est inférieure à celle de l'algorithme naïf, mais il est plus rapide dans les cas où $n>100$ environ et apparaît dans plusieurs bibliothèques logicielles, comme BLAS. Il est très utile pour les grandes matrices sur des domaines exacts tels que les champs finis, où la stabilité numérique n'est pas un problème.

L'algorithme en $O(n^{k})$ actuel avec l'exposant k le plus bas est une généralisation de l'algorithme Coppersmith – Winograd qui a une complexité asymptotique de $O(n^{2,3728639})$, proposé par François Le Gall. L'algorithme de Le Gall et l'algorithme de Coppersmith – Winograd sur lequel il est basé sont similaires à l'algorithme de Strassen : on trouve un moyen de multiplier deux matrices $k\times k$ avec moins de $k^{3}$ multiplications, puis cette technique est appliquée récursivement. Cependant, le coefficient constant caché par la notation grand O est si grand que ces algorithmes ne valent que pour les matrices trop grandes pour être manipulées sur les ordinateurs actuels en pratique.

Étant donné que tout algorithme de multiplication de deux matrices $n\times n$ doit traiter l'ensemble des $2n^{2}$ entrées, il existe une borne inférieure asymptotique des opérations qui est $\Omega (n^{2})$ . Raz a démontré une borne inférieure de $\Omega (n^{2}\log n)$ pour les circuits arithmétiques à coefficient borné sur les nombres réels ou complexes.

Cohn et al. placent les méthodes telles que les algorithmes Strassen et Coppersmith – Winograd dans un contexte de théorie des groupes entièrement différent, en utilisant des triplets de sous-ensembles de groupes finis qui satisfont une propriété de disjonction appelée propriété du triple produit. Ils montrent que si les familles de produits en couronne de groupes abéliens avec des groupes symétriques réalisent des familles de sous-ensembles triples avec une version simultanée du triple produit, alors il existe des algorithmes de multiplication matricielle avec une complexité essentiellement quadratique. La plupart des chercheurs pensent que c'est effectivement le cas. Cependant, Alon, Shpilka et Chris Umans ont récemment montré que certaines de ces conjectures impliquant une multiplication matricielle rapide sont incompatibles avec une autre conjecture plausible, la conjecture de tournesol.

L'algorithme de Freivalds est un algorithme de Monte Carlo simple qui, étant donné les matrices $A$, $B$ et $C$, vérifie en un temps $\Theta (n^{2})$ si $AB=C$.

Algorithmes distribués et évitant la communication

Sur les architectures modernes à mémoire hiérarchique, le coût de chargement et de stockage des éléments de matrice d'entrée a tendance à dominer le coût des calculs arithmétiques. Sur une seule machine, il s'agit de la quantité de données transférées entre la RAM et le cache, tandis que dans un système multi-nœuds à mémoire distribuée, il s'agit plutôt de la quantité de données transférée entre les nœuds. Dans les deux cas, il s'agit de la bande passante. L'algorithme naïf utilisant trois boucles imbriquées utilise une bande passante de communication $\Omega (n^{3})$.

L'algorithme de Cannon, également connu sous le nom d'algorithme 2D, est un algorithme de limitation des mouvements de données qui partitionne chaque matrice d'entrée en une matrice de blocs dont les éléments sont des sous-matrices de taille ${\sqrt {M/3}}\times {\sqrt {M/3}}$, où $M$ est la taille de la mémoire rapide. L'algorithme naïf est ensuite utilisé sur les matrices de blocs, calculant les produits de sous-matrices intégralement en mémoire rapide. Cela réduit la bande passante à $O\left({\frac {n^{3}}{\sqrt {M}}}\right)$, ce qui est asymptotiquement optimal (pour les algorithmes effectuant le calcul de $\Omega (n^{3})$).

Dans un environnement distribué avec $p$ processeurs disposés dans un réseau maillé 2D de forme ${\sqrt {p}}\times {\sqrt {p}}$, une sous-matrice du résultat peut être affectée à chaque processeur, et le produit peut être calculé avec chaque processeur transmettant $O\left({\frac {n^{2}}{\sqrt {p}}}\right)$ mots, ce qui est asymptotiquement optimal si l'on suppose que chaque nœud stocke les O(n²/p) éléments minimaux. Ceci peut être amélioré par l'algorithme 3D, qui organise les processeurs en une maille de cube 3D, affectant chaque produit de deux sous-matrices d'entrée à un seul processeur. Les sous-matrices résultantes sont ensuite générées en effectuant une réduction sur chaque ligne. Cet algorithme transmet $O\left({\frac {n^{2}}{p^{2/3}}}\right)$ mots par processeur, ce qui est asymptotiquement optimal. Cependant, cela nécessite de répliquer chaque élément de matrice d'entrée $p^{1/3}$ fois, et nécessite donc un facteur de $p^{1/3}$ mémoire de plus que ce qui est nécessaire pour stocker les entrées. Cet algorithme peut être combiné avec celui de Strassen pour réduire davantage le temps d'exécution. Les algorithmes dits "2,5D" offrent un compromis continu entre l'utilisation de la mémoire et la bande passante. Sur les environnements informatiques distribués modernes tels que MapReduce, des algorithmes de multiplication spécialisés ont été développés.

Algorithmes pour les maillages

Multiplication matricielle terminée en $2n-1$ étapes pour deux matrices $n\times n$ sur un maillage croisé.

Il existe une variété d'algorithmes de multiplication sur les réseaux maillés. Pour la multiplication de deux matrices de taille $n\times n$ sur un réseau maillé bidimensionnel standard en utilisant l'algorithme de Cannon 2D, on peut effectuer la multiplication en $3n-2$ étapes. Le tableau standard est inefficace car les données des deux matrices n'arrivent pas simultanément et doivent être remplies de zéros.

Le résultat est encore plus rapide sur un maillage croisé à deux couches, où seulement $2n-1$ étapes sont nécessaires. Les performances s'améliorent encore pour les calculs répétés conduisant à une efficacité de 100%.

Voir également

• Algorithme CYK, algorithme de §Valiant
• Multiplication de chaînes matricielles

Références

[[Catégorie:Pages avec des traductions non relues]]