Pourquoi MATLAB est-il si rapide en multiplication matricielle?

Je fais des benchmarks avec CUDA, C ++, C #, Java et j'utilise MATLAB pour la vérification et la génération de matrices. Lorsque j'effectue une multiplication matricielle avec MATLAB, des 2048x2048matrices encore plus grandes sont multipliées presque instantanément.

             1024x1024   2048x2048   4096x4096
             ---------   ---------   ---------
CUDA C (ms)      43.11      391.05     3407.99
C++ (ms)       6137.10    64369.29   551390.93
C# (ms)       10509.00   300684.00  2527250.00
Java (ms)      9149.90    92562.28   838357.94
MATLAB (ms)      75.01      423.10     3133.90

Seul CUDA est compétitif, mais je pensais qu'au moins C ++ serait un peu proche et pas 60 fois plus lent. Je ne sais pas non plus quoi penser des résultats C #. L'algorithme est exactement le même que C ++ et Java, mais il y a un bond 2048de géant 1024.

Comment MATLAB effectue-t-il une multiplication matricielle si rapidement?

Code C ++:

float temp = 0;
timer.start();
for(int j = 0; j < rozmer; j++)
{
    for (int k = 0; k < rozmer; k++)
    {
        temp = 0;
        for (int m = 0; m < rozmer; m++)
        {
            temp = temp + matice1[j][m] * matice2[m][k];
        }
        matice3[j][k] = temp;
    }
}
timer.stop();

Voici mes résultats en utilisant MATLAB R2011a + Parallel Computing Toolbox sur une machine avec une Tesla C2070:

>> A = rand(1024); gA = gpuArray(A);
% warm up by executing the operations a couple of times, and then:
>> tic, C = A * A; toc
Elapsed time is 0.075396 seconds.
>> tic, gC = gA * gA; toc
Elapsed time is 0.008621 seconds.

MATLAB utilise des bibliothèques hautement optimisées pour la multiplication matricielle, raison pour laquelle la multiplication matricielle simple MATLAB est si rapide. La gpuArrayversion utilise MAGMA .

Mise à jour à l'aide de R2014a sur une machine avec un Tesla K20c, et les nouvelles fonctions timeitet gputimeit:

>> A = rand(1024); gA = gpuArray(A);
>> timeit(@()A*A)
ans =
    0.0324
>> gputimeit(@()gA*gA)
ans =
    0.0022

Mise à jour à l'aide de R2018b sur une machine WIN64 avec 16 cœurs physiques et un Tesla V100:

>> timeit(@()A*A)
ans =
    0.0229
>> gputimeit(@()gA*gA)
ans =
   4.8019e-04

(NB: à un moment donné (j'oublie quand exactement) est gpuArraypassé de MAGMA à cuBLAS - MAGMA est encore utilisé pour certaines gpuArrayopérations)

Ce genre de question est récurrent et devrait recevoir une réponse plus claire que "MATLAB utilise des bibliothèques hautement optimisées" ou "MATLAB utilise le MKL" pour une fois sur Stack Overflow.

L'histoire:

La multiplication matricielle (avec matrice-vecteur, multiplication vecteur-vecteur et plusieurs des décompositions matricielles) est (sont) les problèmes les plus importants en algèbre linéaire. Les ingénieurs résolvent ces problèmes avec les ordinateurs depuis les premiers jours.

Je ne suis pas un expert en histoire, mais apparemment à l'époque, tout le monde vient de réécrire sa version FORTRAN avec de simples boucles. Une certaine normalisation est alors venue avec l'identification des «noyaux» (routines de base) dont la plupart des problèmes d'algèbre linéaire avaient besoin pour être résolus. Ces opérations de base ont ensuite été standardisées dans une spécification appelée: Basic Linear Algebra Subprograms (BLAS). Les ingénieurs pourraient alors appeler ces routines BLAS standard et bien testées dans leur code, facilitant ainsi leur travail.

BLAS:

BLAS a évolué du niveau 1 (la première version définissant les opérations vectorielles scalaires et vectorielles) au niveau 2 (opérations vectorielles-matrices) au niveau 3 (opérations matrices-matrices), et a fourni de plus en plus de «noyaux» donc standardisés plus et plus des opérations fondamentales d'algèbre linéaire. Les implémentations originales de FORTRAN 77 sont toujours disponibles sur le site Web de Netlib .

Vers de meilleures performances:

Ainsi, au fil des années (notamment entre les versions BLAS niveau 1 et niveau 2: début des années 80), le matériel a changé, avec l'avènement des opérations vectorielles et des hiérarchies de cache. Ces évolutions ont permis d'augmenter considérablement les performances des sous-programmes BLAS. Différents éditeurs sont alors venus avec leur implémentation de routines BLAS qui étaient de plus en plus efficaces.

Je ne connais pas toutes les implémentations historiques (je ne suis pas né ou enfant à l'époque), mais deux des plus notables sont sorties au début des années 2000: l'Intel MKL et le GotoBLAS. Votre Matlab utilise l'Intel MKL, qui est un très bon BLAS optimisé, et cela explique les excellentes performances que vous voyez.

Détails techniques sur la multiplication matricielle:

Alors pourquoi Matlab (le MKL) est-il si rapide à dgemm(multiplication matrice-matrice générale à double précision)? En termes simples: parce qu'il utilise la vectorisation et une bonne mise en cache des données. En termes plus complexes: voir l' article fourni par Jonathan Moore.

Fondamentalement, lorsque vous effectuez votre multiplication dans le code C ++ que vous avez fourni, vous n'êtes pas du tout compatible avec le cache. Puisque je soupçonne que vous avez créé un tableau de pointeurs vers des tableaux de lignes, vos accès dans votre boucle interne à la k-ème colonne de "matice2": matice2[m][k]sont très lents. En effet, lorsque vous y accédez matice2[0][k], vous devez récupérer le k-ième élément du tableau 0 de votre matrice. Ensuite, dans l'itération suivante, vous devez accéder matice2[1][k], qui est le k-ième élément d'un autre tableau (le tableau 1). Ensuite, dans l'itération suivante, vous accédez à un autre tableau, et ainsi de suite ... Comme la matrice entière matice2ne peut pas tenir dans les caches les plus élevés (ses 8*1024*1024octets sont grands), le programme doit récupérer l'élément souhaité de la mémoire principale, perdant beaucoup de temps.

Si vous venez de transposer la matrice, afin que les accès soient dans des adresses mémoire contiguës, votre code s'exécuterait déjà beaucoup plus rapidement car maintenant le compilateur peut charger des lignes entières dans le cache en même temps. Essayez simplement cette version modifiée:

timer.start();
float temp = 0;
//transpose matice2
for (int p = 0; p < rozmer; p++)
{
    for (int q = 0; q < rozmer; q++)
    {
        tempmat[p][q] = matice2[q][p];
    }
}
for(int j = 0; j < rozmer; j++)
{
    for (int k = 0; k < rozmer; k++)
    {
        temp = 0;
        for (int m = 0; m < rozmer; m++)
        {
            temp = temp + matice1[j][m] * tempmat[k][m];
        }
        matice3[j][k] = temp;
    }
}
timer.stop();

Vous pouvez donc voir à quel point la seule localité du cache a considérablement augmenté les performances de votre code. Désormais, les dgemmimplémentations réelles exploitent cela à un niveau très étendu: elles effectuent la multiplication sur des blocs de la matrice définie par la taille du TLB (Translation lookaside buffer, long story short: what can effective cached), de sorte qu'ils streament vers le processeur exactement la quantité de données qu'il peut traiter. L'autre aspect est la vectorisation, ils utilisent les instructions vectorisées du processeur pour un débit d'instructions optimal, ce que vous ne pouvez pas vraiment faire à partir de votre code C ++ multiplateforme.

Enfin, les personnes prétendant que c'est à cause de l'algorithme de Strassen ou Coppersmith – Winograd ont tort, ces deux algorithmes ne sont pas implémentables en pratique, à cause des considérations matérielles mentionnées ci-dessus.

Voilà pourquoi . MATLAB n'effectue pas de multiplication de matrice naïve en boucle sur chaque élément comme vous l'avez fait dans votre code C ++.

Bien sûr, je suppose que vous venez d'utiliser C=A*Bau lieu d'écrire vous-même une fonction de multiplication.

Matlab a incorporé LAPACK il y a quelque temps, donc je suppose que leur multiplication matricielle utilise quelque chose au moins aussi rapide. Le code source et la documentation de LAPACK sont facilement disponibles.

Vous pouvez également consulter l'article de Goto et Van De Geijn "Anatomy of High-Performance Matrix Multiplication" sur http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.140.1785&rep=rep1&type=pdf

La réponse est que les bibliothèques LAPACK et BLAS rendent MATLAB incroyablement rapide lors des opérations matricielles, pas de code propriétaire par les gens de MATLAB.

Utilisez les bibliothèques LAPACK et / ou BLAS dans votre code C ++ pour les opérations matricielles et vous devriez obtenir des performances similaires à celles de MATLAB. Ces bibliothèques devraient être disponibles gratuitement sur n'importe quel système moderne et certaines parties ont été développées au fil des décennies dans le milieu universitaire. Notez qu'il existe plusieurs implémentations, y compris certaines sources fermées telles que Intel MKL .

Une discussion sur la façon dont BLAS obtient des performances élevées est disponible ici.

BTW, c'est une douleur sérieuse dans mon expérience d'appeler les bibliothèques LAPACK directement à partir de c (mais ça vaut le coup). Vous devez lire la documentation TRÈS précisément.

Lorsque vous effectuez une multiplication matricielle, vous utilisez une méthode de multiplication naïve qui prend du temps O(n^3).

Il existe un algorithme de multiplication matricielle qui prend O(n^2.4). Ce qui signifie que n=2000votre algorithme nécessite environ 100 fois plus de calculs que le meilleur algorithme.
Vous devriez vraiment consulter la page wikipedia pour la multiplication matricielle pour plus d'informations sur les moyens efficaces de l'implémenter.

Selon votre version de Matlab, je pense qu'il utilise déjà votre GPU.

Autre chose; Matlab garde une trace de nombreuses propriétés de votre matrice; que ce soit sa diagonale, hermétienne, etc., et spécialise ses algorithmes basés sur celle-ci. Peut-être que sa spécialisation est basée sur la matrice zéro que vous lui transmettez, ou quelque chose comme ça? Peut-être est-ce la mise en cache des appels de fonction répétés, ce qui gâche votre minutage? Peut-être optimise-t-il les produits matriciels inutilisés répétés?

Pour éviter que de telles choses ne se produisent, utilisez une matrice de nombres aléatoires et assurez-vous de forcer l'exécution en imprimant le résultat sur l'écran ou le disque ou quelque chose de semblable.

MATLAB utilise une implémentation hautement optimisée de LAPACK d'Intel connue sous le nom de Intel Math Kernel Library (Intel MKL) - en particulier la fonction dgemm . La vitesse Cette bibliothèque tire parti des fonctionnalités du processeur, notamment des instructions SIMD et des processeurs multicœurs. Ils ne documentent pas l'algorithme spécifique qu'ils utilisent. Si vous appelez Intel MKL à partir de C ++, vous devriez voir des performances similaires.

Je ne sais pas quelle bibliothèque MATLAB utilise pour la multiplication GPU, mais probablement quelque chose comme nVidia CUBLAS .

La réponse générale à "Pourquoi matlab est-il plus rapide à faire xxx que d'autres programmes" est que matlab a beaucoup de fonctions intégrées et optimisées.

Les autres programmes utilisés n'ont souvent pas ces fonctions, de sorte que les gens appliquent leurs propres solutions créatives, qui sont étonnamment plus lentes que le code optimisé professionnellement.

Cela peut être interprété de deux manières:

1) La méthode commune / théorique: Matlab n'est pas beaucoup plus rapide, vous faites simplement le mauvais test

2) La manière réaliste: Pour ce truc, Matlab est plus rapide en pratique parce que les langages comme le c ++ sont trop facilement utilisés de manière inefficace.

Le contraste net n'est pas seulement dû à l'incroyable optimisation de Matlab (comme déjà discuté par de nombreuses autres réponses), mais aussi à la façon dont vous avez formulé la matrice en tant qu'objet.

Il semble que vous ayez fait une liste de listes à matrice? Une liste de listes contient des pointeurs vers des listes qui contiennent ensuite vos éléments de matrice. Les emplacements des listes contenues sont attribués arbitrairement. Lorsque vous parcourez votre premier index (numéro de ligne?), Le temps d'accès à la mémoire est très important. En comparaison, pourquoi n'essayez-vous pas d'implémenter la matrice sous forme de liste / vecteur unique en utilisant la méthode suivante?

#include <vector>

struct matrix {
    matrix(int x, int y) : n_row(x), n_col(y), M(x * y) {}
    int n_row;
    int n_col;
    std::vector<double> M;
    double &operator()(int i, int j);
};

Et

double &matrix::operator()(int i, int j) {
    return M[n_col * i + j];
}

Le même algorithme de multiplication doit être utilisé pour que le nombre de flop soit le même. (n ^ 3 pour les matrices carrées de taille n)

Je vous demande de le chronométrer pour que le résultat soit comparable à ce que vous aviez auparavant (sur la même machine). Avec la comparaison, vous montrerez exactement à quel point le temps d'accès à la mémoire peut être important!

C'est lent en C ++ car vous n'utilisez pas le multithreading. Essentiellement, si A = BC, où ce sont toutes des matrices, la première ligne de A peut être calculée indépendamment de la 2ème ligne, etc. Si A, B et C sont toutes des matrices n par n, vous pouvez accélérer la multiplication de un facteur de n ^ 2, comme

a_ {i, j} = somme_ {k} b_ {i, k} c_ {k, j}

Si vous utilisez, par exemple, Eigen [ http://eigen.tuxfamily.org/dox/GettingStarted.html ], le multithreading est intégré et le nombre de threads est réglable.

Parce que MATLAB est un langage de programmation initialement développé pour l'algèbre linéaire numérique (manipulations matricielles), qui possède des bibliothèques spécialement développées pour les multiplications matricielles. Et maintenant, MATLAB peut également utiliser les GPU (Graphics Processing Unit) pour cela en plus.

Et si nous regardons vos résultats de calcul:

             1024x1024   2048x2048   4096x4096
             ---------   ---------   ---------
CUDA C (ms)      43.11      391.05     3407.99
C++ (ms)       6137.10    64369.29   551390.93
C# (ms)       10509.00   300684.00  2527250.00
Java (ms)      9149.90    92562.28   838357.94
MATLAB (ms)      75.01      423.10     3133.90

alors nous pouvons voir que non seulement MATLAB est si rapide en multiplication matricielle: CUDA C (langage de programmation de NVIDIA) a de meilleurs résultats que MATLAB. CUDA C a également des bibliothèques spécialement développées pour les multiplications matricielles et utilise les GPU.

Brève histoire de MATLAB

Cleve Moler, président du département d'informatique de l'Université du Nouveau-Mexique, a commencé à développer MATLAB à la fin des années 1970. Il l'a conçu pour permettre à ses étudiants d'accéder à LINPACK (une bibliothèque de logiciels pour effectuer l'algèbre linéaire numérique) et EISPACK(est une bibliothèque de logiciels pour le calcul numérique de l'algèbre linéaire) sans qu'ils aient à apprendre Fortran. Il s'est rapidement répandu dans d'autres universités et a trouvé un public important au sein de la communauté des mathématiques appliquées. Jack Little, un ingénieur, y a été exposé lors d'une visite de Moler à l'Université de Stanford en 1983. Reconnaissant son potentiel commercial, il s'est joint à Moler et Steve Bangert. Ils ont réécrit MATLAB en C et ont fondé MathWorks en 1984 pour poursuivre son développement. Ces bibliothèques réécrites étaient connues sous le nom de JACKPAC. En 2000, MATLAB a été réécrit pour utiliser un nouvel ensemble de bibliothèques pour la manipulation de matrices, LAPACK (est une bibliothèque logicielle standard pour l'algèbre linéaire numérique).

La source

Qu'est-ce que CUDA C

CUDA C utilise également des bibliothèques spécialement développées pour les multiplications matricielles comme OpenGL (Open Graphics Library). Il utilise également GPU et Direct3D (sur MS Windows).

La plate-forme CUDA est conçue pour fonctionner avec des langages de programmation tels que C, C ++ et Fortran. Cette accessibilité permet aux spécialistes de la programmation parallèle d'utiliser plus facilement les ressources GPU, contrairement aux API antérieures comme Direct3D et OpenGL , qui nécessitaient des compétences avancées en programmation graphique. En outre, CUDA prend en charge les cadres de programmation tels que OpenACC et OpenCL .

Exemple de flux de traitement CUDA:

1. Copier les données de la mémoire principale vers la mémoire GPU
2. Le processeur lance le noyau de calcul GPU
3. Les cœurs CUDA du GPU exécutent le noyau en parallèle
4. Copiez les données résultantes de la mémoire GPU vers la mémoire principale

Comparaison des vitesses d'exécution du processeur et du GPU

Nous avons exécuté un benchmark dans lequel nous avons mesuré le temps nécessaire pour exécuter 50 pas de temps pour des tailles de grille de 64, 128, 512, 1024 et 2048 sur un processeur Intel Xeon X5650, puis en utilisant un GPU NVIDIA Tesla C2050.

Pour une taille de grille de 2048, l'algorithme montre une diminution de 7,5 fois du temps de calcul de plus d'une minute sur le CPU à moins de 10 secondes sur le GPU. Le graphique à l'échelle logarithmique montre que le processeur est en fait plus rapide pour les petites tailles de grille. Au fur et à mesure que la technologie évolue et mûrit, les solutions GPU sont de plus en plus capables de gérer des problèmes plus petits, une tendance que nous prévoyons de continuer.

La source

De l'introduction au guide de programmation CUDA C:

Poussé par la demande insatiable du marché pour des graphiques 3D haute définition en temps réel, l'unité de processeur graphique programmable ou GPU a évolué pour devenir un processeur multithread hautement parallèle avec une puissance de calcul énorme et une bande passante mémoire très élevée, comme illustré par Figure 1et Figure 2.

Figure 1. Opérations en virgule flottante par seconde pour le CPU et le GPU

Graphique 2 . Bande passante mémoire pour le CPU et le GPU

La raison de l'écart de capacité en virgule flottante entre le CPU et le GPU est que le GPU est spécialisé pour le calcul intensif en calcul et hautement parallèle - exactement ce qu'est le rendu graphique - et donc conçu de manière à ce que davantage de transistors soient consacrés au traitement des données. plutôt que la mise en cache des données et le contrôle de flux, comme illustré schématiquement par Figure 3.

Graphique 3 . Le GPU consacre plus de transistors au traitement des données

Plus précisément, le GPU est particulièrement bien adapté pour résoudre les problèmes qui peuvent être exprimés sous forme de calculs parallèles aux données - le même programme est exécuté sur de nombreux éléments de données en parallèle - avec une intensité arithmétique élevée - le rapport entre les opérations arithmétiques et les opérations de mémoire. Parce que le même programme est exécuté pour chaque élément de données, il y a moins d'exigences pour un contrôle de flux sophistiqué, et comme il est exécuté sur de nombreux éléments de données et a une intensité arithmétique élevée, la latence d'accès à la mémoire peut être masquée avec des calculs au lieu de caches de Big Data .

Le traitement parallèle aux données mappe les éléments de données sur des threads de traitement parallèles. De nombreuses applications qui traitent de grands ensembles de données peuvent utiliser un modèle de programmation parallèle aux données pour accélérer les calculs. Dans le rendu 3D, de grands ensembles de pixels et de sommets sont mappés sur des threads parallèles. De même, les applications de traitement d'images et de médias telles que le post-traitement des images rendues, le codage et le décodage vidéo, la mise à l'échelle d'image, la vision stéréo et la reconnaissance de formes peuvent mapper des blocs d'image et des pixels sur des fils de traitement parallèles. En fait, de nombreux algorithmes en dehors du domaine du rendu et du traitement d'images sont accélérés par le traitement parallèle aux données, du traitement général du signal ou de la simulation physique à la finance informatique ou à la biologie computationnelle.

La source

Lecture avancée

• GPU (unité de traitement graphique)
• MATLAB
• Guide de programmation CUDA C
• Utilisation de GPU dans MATLAB

• Sous-programmes d'algèbre linéaire de base (BLAS)

• Anatomy of High-Performance Matrix Multiplication , de Kazushige Goto et Robert A. Van De Geijn

Quelques facs intéressants

J'ai écrit une multiplication matricielle C ++ qui est aussi rapide que celle de Matlab mais cela a pris quelques précautions. (Avant que Matlab utilisait des GPU pour cela).

Сitation de cette réponse .