Surpasser la multiplication de matrices de NumPy avec 150 lignes de C
(salykova.github.io)- La multiplication de matrices de NumPy s’appuie sur des bibliothèques BLAS externes, mais cette implémentation vise à porter les performances mono-thread et multi-thread au niveau de BLAS avec seulement du C pur et FMA3·AVX2
- Le cœur des performances consiste à découper $C$ en petits blocs et à utiliser un microkernel 16×6 qui répète des rank-1 updates dans les registres
YMMafin de réduire les accès mémoire - Pour les matrices de taille arbitraire, la gestion des bords peut vite devenir un goulot d’étranglement ; l’implémentation combine donc des stores masqués et des tampons avec padding à 0 pour éviter la baisse de performance des loads masqués
- La réutilisation du cache est assurée par un blocking
k_c,m_c,n_c, et les performances maximales réelles dépendent fortement du tuning du nombre de threads, de la taille du kernel et de la taille des tuiles - AVX-512 a été exclu afin de prendre en charge un éventail plus large de CPU ; sur les CPU AVX-512, BLAS peut donc être plus rapide, et la comparaison avec OpenBLAS a elle aussi été réalisée avec AVX-512 désactivé
Objectif de l’implémentation et base de comparaison
- Le code de l’implémentation est publié dans sgemm.c et optimise la multiplication de matrices FP32 multi-thread sur les processeurs récents
- NumPy dépend de bibliothèques BLAS externes pour les opérations d’algèbre linéaire comme la multiplication de matrices
- Exemples : Intel MKL, Accelerate, BLIS, GotoBLAS, OpenBLAS
- OpenBLAS, GotoBLAS et BLIS sont écrits en C/FORTRAN/Assembly et incluent des implémentations de multiplication de matrices optimisées manuellement pour chaque microarchitecture CPU
- L’objectif est une implémentation de multiplication de matrices écrite en C pur, sans assembleur bas niveau, tout en respectant les conditions suivantes
- Fonctionner avec des matrices de tailles arbitraires
- S’exécuter sur les processeurs x86-64 récents
- Rivaliser avec les bibliothèques BLAS existantes
- Rester simple et facile à étendre
- Les références sont Fast Multidimensional Matrix Multiplication on CPU from Scratch de Simon Boehm, Matrix Multiplication de Sergey Slotin, Can you multiply a matrix? de Geohot, ainsi que des articles sur GotoBLAS et BLIS
Conditions de benchmark et calcul des FLOPS
- L’environnement de test est AMD Ryzen 7 9700X, 32 Go DDR5 6000 MHz CL36, OpenBLAS 0.3.26, GCC 13.3, Ubuntu 24.04.1 LTS
- Les flags de compilation utilisés sont
-O3 -march=native -mno-avx512f -fopenmp - Pour une comparaison équitable, il faut définir le
TARGETapproprié lors de l’installation d’OpenBLAS et désactiver les instructions AVX-512- Les processeurs Zen4/5 se compilent avec
make TARGET=ZEN - Sinon, OpenBLAS utilise les instructions AVX-512 par défaut
- Les processeurs Zen4/5 se compilent avec
- La multiplication de matrices FP32 d’OpenBLAS est exécutée via l’API
cblas_sgemm - Le benchmark porte sur des matrices carrées
- Évaluation de
m=n=k=200àm=n=k=10000, par pas de200 - La multiplication de matrices est répétée
n_iterfois, et le temps d’exécution médian est utilisé pour mesurer les performances
- Évaluation de
- Multiplier une matrice $A$ de taille $M \times K$ par une matrice $B$ de taille $K \times N$ représente un total de $2MNK$ FLOP
- Les performances sont calculées avec
FLOPS=(2*m*n*k)/exec_time
- Les performances sont calculées avec
Limites théoriques et base SIMD
- Les CPU x86-64 récents traitent plusieurs données en parallèle grâce aux extensions SIMD
- Les principales instructions sont AVX2 et FMA
- Toutes deux utilisent des registres
YMMde 256 bits - Chaque registre
YMMpeut contenir 8 floats de 32 bits
- Toutes deux utilisent des registres
- L’instruction FMA
VFMADD231PSeffectue une opération packed single de la formeYMM1 = YMM2 * YMM3 + YMM1 - Sur Ryzen 9700X, le débit du fused multiply-add est de 0,5 cycle/instruction, soit 2 instructions par cycle
- En théorie, le Ryzen 9700X peut effectuer 32 FLOP par cycle sur un seul cœur
- Le calcul est
8 floats × 2(add+mul) × 2(1/TP) - En supposant une fréquence soutenue de 4,7 GHz sur 8 cœurs, le pic théorique multi-thread est estimé à 1203 FLOPS
- Le calcul est
Implémentation de base et microkernel
- Les matrices sont stockées en ordre column-major
A[row][col]est accessible, depuis un pointeur C, viaptr[col*M + row]
- L’implémentation la plus simple parcourt toutes les lignes et colonnes de $C$ et calcule, pour chaque élément, le produit scalaire d’une ligne de $A$ et d’une colonne de $B$
- Le cœur d’une implémentation haute performance est le microkernel, qui découpe $C$ en sous-matrices $m_R \times n_R$ et calcule efficacement chacune d’elles
- Le kernel initialise $\bar{C}$ à zéro dans les registres, puis itère le long de la dimension $K$
- Il charge dans les registres le vecteur colonne de $\bar{A}$ et le vecteur ligne de $\bar{B}$
- Il calcule le produit extérieur de ces deux vecteurs et l’ajoute aux accumulateurs $\bar{C}$
- Chaque étape est une rank-1 update
- Par rapport aux $2K m_R n_R$ accès mémoire de l’approche naïve, cette méthode réduit à $(m_R+n_R)K$ le nombre d’éléments chargés dans les registres
- Les CPU AVX disposent de 16 registres YMM, la taille du kernel doit donc respecter la contrainte suivante
- $(m_R/8) \cdot n_R + m_R/8 + 1 \le 16$
- $m_R$ doit être un multiple de 8
- En théorie, plus $m_R$ et $n_R$ sont grands et proches l’un de l’autre, plus la réduction des accès mémoire est importante, mais en pratique, sur Ryzen 9700X, le kernel 16×6 donne les meilleures performances
- L’implémentation utilise les intrinsics de
immintrin.h__m256est un type vectoriel 256 bits représentant le contenu d’un registreYMM_mm256_loadu_pscharge un vecteur colonne deA_mm256_broadcast_ssdiffuse une valeur scalaire deBdans un vecteur de 8 floats_mm256_fmadd_psmet à jour les accumulateurs_mm256_storeu_psstocke le résultat en mémoire
- L’assembleur généré contient des instructions SIMD FMA comme
vfmadd231psetvbroadcastss
Padding pour matrices de taille arbitraire
- Le kernel 16×6 de base fonctionne directement lorsque $M$ et $N$ sont respectivement des multiples de 16 et de 6
- Sur les zones de bord, si le nombre de colonnes $n$ est inférieur à 6, la boucle de stockage ne s’exécute que jusqu’à
j < n - Lorsque le nombre de lignes $m$ est inférieur à 16,
_mm256_storeu_psstocke 8 éléments à la fois ; un store masqué est donc nécessaire_mm256_maskstore_psn’écrit en mémoire que les éléments dont le bit de masque est activé- Le masque est généré selon le nombre de lignes recouvrées $m$
- Traiter aussi les loads de bord avec
_mm256_maskload_pspeut fortement dégrader les performances du kernel- Les instructions supplémentaires de calcul du masque ajoutent un surcoût
- Comme $n$ n’est pas une constante connue à la compilation, le compilateur a du mal à dérouler efficacement la boucle
- À la place, si $m \neq m_R$, $\bar{A}$ est copié dans un tampon et complété par des zéros ; si $n \neq n_R$, $\bar{B}$ est également copié dans un tampon et rempli de zéros
- L’implémentation correspondante se trouve dans matmul_pad.h
Blocking de cache et réutilisation des données
- Entre les registres et la DRAM se trouve la hiérarchie de cache CPU ; les CPU de bureau récents utilisent généralement des caches L1, L2 et L3
- Le cache est plus rapide que la DRAM mais sa capacité est limitée, il est donc impossible de placer l’intégralité de $A$, $B$ et $C$ en cache
- La méthode consistant à découper les matrices en petits blocs, à les charger dans le cache et à réutiliser les mêmes données pour plusieurs rank-1 updates s’appelle blocking de cache ou tiling
- Le blocking de cache mono-thread prend la forme d’une boucle quintuple similaire à la structure de BLIS
- La boucle la plus externe parcourt la dimension $N$ pour créer les blocs $C_j$ et $B_j$
- La boucle suivante parcourt la dimension $K$ pour créer les blocs $A_j$ et $B_p$
- $B_p$ est packé pour devenir $\tilde{B}_p$ et, si nécessaire, complété par des zéros afin de favoriser sa réutilisation dans le cache L3
- La boucle suivante parcourt la dimension $M$ pour créer les blocs $C_i$ et $A_j$, et $A_j$ est packé pour devenir $\tilde{A}_j$
- Les deux dernières boucles découpent les blocs de cache en panneaux $m_R \times k_c$ et $k_c \times n_R$, puis les transmettent au kernel
- Les matrices packées $\tilde{A}_j$ et $\tilde{B}_p$ sont stockées différemment
- Les panneaux internes de $\tilde{A}_j$ sont stockés en column-major
- Les panneaux internes de $\tilde{B}_p$ sont stockés en row-major
- Les paramètres de blocking de cache doivent être ajustés aux tailles de cache de chaque modèle de CPU
- $k_c \times n_c$ est un point de départ pour remplir le cache L3
- $m_c \times k_c$ est un point de départ pour remplir le cache L2
- $k_c \times n_R$ est un point de départ pour remplir le cache L1
- En pratique, des valeurs supérieures aux valeurs théoriques donnent souvent de meilleures performances, et comme le CPU gère automatiquement le placement en cache, il faut concevoir les boucles et les schémas d’accès au niveau algorithmique
- L’implémentation se trouve dans matmul_cache.h
Micro-optimisations du kernel
- Au lieu de définir les accumulateurs sous forme de tableau, comme
__m256 C_buffer[6][2], les variables d’accumulateurs sont explicitement déroulées et déclarées séparément - Cette approche aide GCC à mieux optimiser le code et à éviter le register spilling
- Le calcul du masque est également modifié pour utiliser des instructions vectorielles
- Un tableau statique
mask[32]est utilisé avec_mm256_cvtepi8_epi32et_mm_loadu_si64
- Un tableau statique
- Cette implémentation se trouve dans matmul_micro.h
Stratégie multi-threading
- La parallélisation vise à la fois les opérations arithmétiques et le packing
- Les 5e, 4e et 3e boucles en partant du microkernel répètent des blocs de la taille des blocs de cache
- Pour occuper tous les threads, le nombre d’itérations doit être au moins égal au nombre de threads
- Les dimensions des matrices d’entrée doivent être approximativement au moins égales à
nombre de threads × taille de bloc de cache
- Sur Ryzen 9700X, les tailles de blocs de cache qui donnent de bonnes performances en mono-thread sont $n_c=1535$, $m_c=1024$
- Pour exploiter les 8 cœurs, il faut au minimum une dimension de taille $\max(m_c,n_c) \times 8 = 1535 \times 8 = 12280$
- À l’inverse, les deux dernières boucles répètent de petits blocs $m_R$, $n_R$ et conviennent donc à la parallélisation
- En général, $m_R$ et $n_R$ sont inférieurs à 20
- Choisir $m_c$ et $n_c$ comme multiples du nombre de cœurs permet de répartir équitablement le travail
- Sur Ryzen 9700X, la meilleure performance est obtenue en parallélisant ensemble les deux boucles internes avec
#pragma omp parallel for collapse(2) num_threads(NTHREADS) - Sur les processeurs avec beaucoup de cœurs, en particulier au-delà de 16 cœurs, on peut envisager le parallélisme imbriqué et la parallélisation de 2 à 3 boucles
- Le packing de $\tilde{A}$ et $\tilde{B}$ est également parallélisé avec OpenMP
pack_blockAest parallélisé en parcourantmcpar pas deMRpack_blockBest parallélisé en parcourantncpar pas deNR
- Les paramètres qui donnent de bonnes performances sur Ryzen 9700X dans l’implémentation multi-thread sont les suivants
- $m_c = m_R \times \text{nombre de threads} \times 5$
- $n_c = n_R \times \text{nombre de threads} \times 50$
- L’implémentation multi-thread finale se trouve dans matmul_parallel.h
1 commentaires
Commentaires sur Hacker News
Si l’idée de cet article est simplement qu’il reste généralement de la marge de performance, alors il sous-estime plutôt l’ampleur possible des améliorations. Et ce, même si l’effort investi dans les bibliothèques de multiplication de matrices est bien plus important que dans la plupart des logiciels.
À moins d’avoir déjà un code fortement optimisé, il est courant d’obtenir, sans énorme effort, des gains de 10 à 1000 fois, voire plus par rapport au code existant. En gros, par ordre d’importance, le plus crucial est de savoir si le choix de l’algorithme est approprié et si l’on peut supprimer purement et simplement le travail à effectuer ; réduire les opérations coûteuses comme les allers-retours vers le noyau ou les
malloccompte aussi beaucoup.Pour la vectorisation, les intrinsics vectoriels explicites sont utiles, mais le simple fait de réorganiser les données, par exemple en passant d’un tableau de structures à une structure de tableaux, produit souvent le même code machine. L’efficacité du cache est également importante, et dans le code parallèle cela devient plus complexe lorsque les données ne sont pas bien isolées par thread, comme avec le faux partage. Enfin, on peut aussi faire des optimisations spécifiques au matériel, via des intrinsics ou de l’assembleur écrit à la main.
En modifiant la requête pour que la jointure se fasse sur le serveur distant et que seules 5 à 10 lignes soient envoyées sur le réseau, c’est devenu immédiatement plus rapide. Il y a toujours un surcoût fixe et de la latence, mais si l’on envoie beaucoup plus de données que nécessaire sur une connexion réseau, les performances finissent par s’effondrer. L’article « It's the latency, stupid », sur l’impact de la latence, vaut aussi la lecture : http://www.stuartcheshire.org/rants/latency.html
Dans l’ensemble, je suis d’accord avec les considérations ci-dessus et leur ordre approximatif.
Beaucoup d’entretiens sont devenus des quiz de mémorisation d’algorithmes obscurs façon « Google fait comme ça », au lieu d’évaluer la capacité à raisonner sur les causes de lenteur d’une implémentation, à la benchmarker et à la corriger.
Les schémas de code courants ne sont pas assez spécialisés pour le matériel et laissent beaucoup de performance sur la table. Cet article en est un exemple intéressant ; une autre démonstration classique est « There's plenty of room at the top ».
https://www.science.org/doi/10.1126/science.aam9744
Pour comprendre ce sujet, les articles du dépôt BLIS sont ce qui se rapproche le plus d’une référence. Je ne vois pas pourquoi on penserait qu’un BLAS optimisé ne donne pas de bonnes performances ; pour des matrices suffisamment grandes, il faut s’attendre à plus de 90 % du pic du CPU.
La dernière fois que j’ai regardé, OpenBLAS en série était globalement comparable à MKL, et BLAS implémente GEMM comme brique de base d’algèbre linéaire, pas
matmul. Je ne comprends pas non plus pourquoi on utiliserait numpy au lieu d’un framework de benchmark, et sur Zen il faudrait à mon avis comparer avec le BLAS d’AMD, c’est-à-dire l’implémentation basée sur BLIS. BLIS avait autrefois une meilleure histoire côté parallélisation qu’OpenBLAS, et AMD BLIS dispose aussi d’un basculement d’implémentation pour les « petites » dimensions ; je ne sais pas si OpenBLAS l’a aujourd’hui.Les intrinsics SIMD ne sont pas indispensables pour vectoriser un micro-kernel ; un bon compilateur C peut vectoriser complètement et dérouler les boucles. Le micro-kernel en C pur de BLIS atteint plus de 80 % des performances de l’implémentation optimisée à la main pour Haswell, avec des tailles de blocs appropriées. La différence vient probablement du préchargement, mais je ne la comprends pas précisément.
Même les compilateurs modernes ne réussissent pas l’auto-vectorisation et le déroulage de boucles correctement dans 100 % des cas.
L’article et l’implémentation ont l’air bons, mais je me demande quel est le « secret ». OpenBLAS est optimisé en assembleur+C depuis des décennies précisément pour ce type de problème ; comment peut-il être battu ?
L’article détaille le cache, etc., mais je me demande si BLAS n’exploite pas ce genre de choses, ou si c’est simplement mieux adapté à un processeur particulier.
L’écart de performance était bien plus marqué sur le débit de pointe que sur le débit moyen, or il y a très peu d’applications où le pic compte vraiment. Le code de benchmark affiché semble faire passer numpy par l’allocateur Python, alors que l’implémentation C ne passe pas par l’allocateur ; c’est donc le premier endroit où chercher une erreur ou une incohérence de microbenchmark. Beaucoup de routines numpy prennent en charge les opérations en place, donc il faudrait explicitement voir un benchmark des versions en place des deux côtés.
numpy effectue aussi des vérifications de bornes et de la gestion d’erreurs indépendantes de l’implémentation sous-jacente, ce qui explique pourquoi il peut être très lent sur de petites matrices, parfois même plus lent que des listes Python brutes. Ajouter quelques milliers de cycles de pur overhead rend vite les choses difficiles à accélérer.
Cette implémentation suit une approche assez méthodique pour saturer les caches concernés et, en un sens, c’est évident, mais des améliorations d’ingénierie claires comme celle-ci méritent d’être mises en avant dans ce type de discussion. OpenBLAS a aussi mobilisé beaucoup de monde, mais il est peu probable qu’ils aient pensé à tout. Pour l’expliquer correctement, il faudrait une analyse approfondie du code des deux côtés.
Pour des implémentations C++ et C, voir les approches par métaprogrammation [2], [3]. Ce qui est vraiment surprenant, c’est que beaucoup de langages modernes comme Matlab, Julia ou Mojo dépendent encore d’OpenBLAS, même s’ils ont bien sûr chacun leurs raisons.
[1] Numeric age for D: Mir GLAS is faster than OpenBLAS and Eigen (2016):
http://blog.mir.dlang.io/glas/benchmark/openblas/2016/09/23/...
[2] Vastly outperforming LAPACK with C++ metaprogramming (2018):
https://wordsandbuttons.online/vastly_outperforming_lapack_w...
[3] Outperforming LAPACK with C metaprogramming (2018):
https://wordsandbuttons.online/outperforming_lapack_with_c_m...
-march=nativepeut donner un avantage, car il compile pour le modèle exact de CPU. numpy a de fortes chances d’avoir été compilé pour une cible x86-64 plus générique et plus ancienne.Sur les CPU Ryzen,
-march=nativeutilise probablement v4, tandis que numpy cible sans doute v1 ou v2.https://en.wikipedia.org/wiki/X86-64#Microarchitecture_level...
L’article est bon, et c’est aussi très bien d’avoir rendu les benchmarks faciles à reproduire. Sur mon Xeon W-2245 16 cœurs à 3,90 GHz,
matmul.ceffectue une multiplication de matrices 8192x8192 en 1,41 s avecgcc -O3et en 1,47 s avecclang -O2, tandis que NumPy prend 1,07 s.Avec un kernel AVX-512, je pense que ce serait beaucoup plus rapide. Une autre raison possible des performances décevantes est OpenMP ; d’après mon expérience, gérer explicitement un pool de threads avec
pthreadspermet de réduire l’overhead. Il vaudrait aussi mieux utilisersysconf(_SC_NPROCESSORS_ONLN)plutôt que de coder en dur le nombre de CPU.Il n’y a pas de raison d’imposer des charges différentes, avec Python d’un côté et C de l’autre. On aurait pu écrire les deux en C, l’un appelant une bibliothèque BLAS et l’autre cette implémentation, pour faire une comparaison à armes égales.
L’overhead n’est pas énorme, mais comme dit ailleurs dans ce fil, il est important de l’appeler correctement. Opposer du code numpy naïf à du code C optimisé n’est clairement pas une comparaison équitable.
Ce n’est pas un chemin critique, mais l’inefficacité de la génération du masque, c’est-à-dire l’utilisation de
bit_mask, me gêne. Des méthodes plus efficaces consisteraient à créer un tableau constant global de la forme{-1,-1,...,0,0,...}et à charger à partir des offsets d’éléments16-m,8-m, ou à comparer un vecteur constant{0,1,2,3,4,...}avecmetm-8diffusés par broadcastCela dit, cela ne concerne qu’une seule colonne de la matrice, et la boucle
maskload/maskstorequi suit prend beaucoup plus de temps, donc c’est vraiment un détail. En particulier, les stores restent lents même sur Zen 4[1], et les instructions AVX-512 sont 6 fois plus rapides alors que la seule différence est qu’elles prennent le masque depuis un registre de masque. Clang vectorise de toute façon automatiquement les décalages, donc ce sera probablement seulement 2 à 3 fois plus lent que ma suggestion[1]: https://uops.info/table.html?search=vmaskmovps&cb_lat=on&cb_...
J’apprécie vraiment les retours qui apportent un nouveau point de vue. Si je me souviens bien, « créer un tableau global de constantes et le charger » était un peu plus lent que le décalage du masque de bits quand je l’avais testé, mais je vais retester pour en être sûr. La méthode consistant à « comparer un vecteur constant
{0, 1, 2, 3, 4, ...}avecmetm-8diffusés par broadcast » est une bonne idée, je vais l’essayerint8_tet, lors du chargement, faire une extension de signe des octets enint32_t. La combinaison_mm_loadu_si64/_mm256_cvtepi8_epi32devrait être compilée en une seule instructionvpmovsxbdavec opérande mémoireAinsi, avec un alignement correct via
alignas(32), tout le tableau de constantes tient dans une seule ligne de cache. Dans le cas d’usage de l’article, il faut deux masques, donc la deuxième instructionvpmovsxbdsera clairement un hit dans le cache L1D, ce qui convient bienEt tinyBLAS de jart ?
https://hacks.mozilla.org/2024/04/llamafiles-progress-four-m...
Et https://justine.lol/matmul/
À part les benchmarks, quelle est la raison de multithreader la multiplication matricielle elle-même ? En pratique, ne vaudrait-il pas mieux utiliser plusieurs threads au niveau de l’algorithme qui utilise la multiplication ?
Mais dans le code HPC, GEMM n’est généralement pas le goulot d’étranglement
Je n’ai encore fait que survoler, mais cet article contient beaucoup de détails et d’explications. Il semble assez bien expliquer comment une multiplication matricielle rapide est implémentée en tenant compte des considérations d’architecture, donc je l’ai ajouté à ma liste de lecture