Pourquoi MATLAB est-il si rapide dans la multiplication matricielle?

ce genre de question est récurrente 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 Stackoverflow.

l'Histoire:

(ainsi que la multiplication matricielle-vectorielle, vectorielle-vectorielle et bon nombre des décompositions matricielles) est (sont) le problème le plus important dans l'algrebra linéaire. Les ingénieurs ont résolu ces problèmes avec les ordinateurs depuis le début.

Je ne suis pas un expert en histoire, mais apparemment à l'époque, tout le monde a réécrit sa version Fortran avec des boucles simples. Une certaine standardisation est alors venue, avec l'identification des" noyaux " (routines de base) que la plupart des problèmes d'algèbre linéaire nécessaire pour être résolu. Ces opérations de base ont ensuite été normalisées dans une spécification appelée: sous-Programmes D'algèbre linéaire de base (BLAS). Les ingénieurs pourraient alors appeler ces, des routines BLAS bien testées dans leur code, rendant leur travail beaucoup plus facile.

BLAS:

BLAS a évolué du niveau 1 (la première version qui a défini les opérations scalaire-vecteur et vecteur-Vecteur) au niveau 2 (opérations vectorielles-matrice) au niveau 3 (opérations matrice-matrice), et a fourni de plus en plus de "noyaux" afin de standardiser de plus en plus les opérations algébriques linéaires fondamentales. Les implémentations originales de Fortran 77 sont: toujours disponible sur Netlib du site internet .

vers une meilleure performance:

ainsi au fil des années (notamment entre les versions BLAS de niveau 1 et de 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 sensiblement les performances des sous-programmes BLAS. Différents vendeurs sont alors venus avec leur mise en œuvre de routines BLAS de plus en plus efficaces.

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

détails techniques sur la multiplication matricielle:

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

en gros, lorsque vous effectuez votre multiplication dans le code C++ que vous avez fourni, vous n'êtes pas du tout compatible avec le cache. Depuis je soupçonne que vous avez créé un tableau de pointeurs rangées de tableaux, vos accès dans votre boucle intérieure à la colonne k-th de "matice2": matice2[m][k] sont très lents. En effet, lorsque vous accédez à matice2[0][k] , vous devez obtenir l'élément k-th du tableau 0 de votre matrice. Puis dans la prochaine itération, vous devez accéder à matice2[1][k] , qui est le k-ième élément d'un autre tableau (tableau 1). Ensuite, dans la prochaine itération, vous accédez encore à un autre tableau, et ainsi de suite... Depuis la matrice entière matice2 ne peut pas tenir dans les caches les plus élevés (il est 8*1024*1024 octets de grande), le programme doit récupérer l'élément désiré de la mémoire principale, perdant beaucoup de temps.

si vous venez de transposer la matrice, de sorte que les accès seraient dans des adresses mémoire contiguës, votre code serait déjà exécuté beaucoup plus rapidement parce que maintenant le compilateur peut charger des lignes entières dans le cache en même temps. Il suffit d'essayer 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();

donc vous pouvez voir comment la localisation de cache a augmenté considérablement les performances de votre code. Maintenant réel Les implémentations dgemm exploitent cela à un niveau très étendu: elles effectuent la multiplication sur des blocs de la matrice définis par la taille du TLB (tampon de translation lookaside, long story short: what can effectively be cached), de sorte qu'elles transmettent au 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'instruction optimal, ce que vous ne pouvez pas vraiment faire à partir de votre code C++ multiplateformes.

enfin, les gens affirmant que c'est à cause de L'algorithme de Strassen ou de Coppersmith–Winograd sont faux, ces deux algorithmes ne sont pas implémentables dans la pratique, en raison de considérations matérielles mentionnées ci-dessus.

C'est pourquoi . MATLAB n'effectue pas une multiplication matricielle naïve en bouclant chaque élément comme vous l'avez fait dans votre code C++.

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

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

vous pouvez aussi regarder L'article de Goto et Van de Geijn "Anatomy of High-Performance Matrix Multiplication" à http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.140.1785&rep=rep1&type=pdf

en faisant la multiplication de matrice, vous utilisez la méthode de multiplication naïve qui prend le temps de O(n^3) .

il existe un algorithme de multiplication matricielle qui prend O(n^2.4) . Ce qui signifie qu'à n=2000 votre algorithme nécessite environ 100 fois plus de calcul que le meilleur algorithme.

Vous devriez vraiment vérifier la page wikipedia pour la multiplication matricielle pour plus d'informations sur les façons efficaces de l'implémenter.

vous devez faire attention à ne pas faire de comparaisons équitables avec C++. Pouvez-vous afficher le code C++ qui montre les boucles internes du noyau que vous utilisez pour la multiplication matricielle? La plupart du temps, je suis préoccupé par votre disposition de la mémoire et si vous faites des choses gaspillées.

j'ai écrit la multiplication matricielle C++ qui est aussi rapide que celle de Matlab mais qui a pris un certain soin. (EDIT: avant que Matlab utilise GPUs pour cela.)

Vous pouvez être pratiquement garanti que Matlab gaspille très peu de cycles sur ces fonctions "intégrées". Ma question Est, où gaspillez-vous vos cycles? (N'en déplaise)

, La réponse est LAPACK et BLAS bibliothèques de MATLAB aveuglément rapide à des opérations matricielles, pas n'importe quel code de propriété 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 librement disponibles sur n'importe quel système moderne et les pièces ont été développées depuis des décennies dans les universités. Notez qu'il y a plusieurs implémentations, y compris certaines sources fermées telles que Intel MKL .

Une discussion de la façon dont BLAS obtient de haute performance , est disponible ici.

BTW, c'est une grave douleur dans mon expérience, à l'appel LAPACK bibliothèques directement à partir de c (mais à peine). Vous devez lire la documentation très précisément.

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

autre chose; Matlab garde trace de beaucoup de propriétés de votre matrice; wether sa diagonale, hermétique, et ainsi de suite, et se spécialise ses algorithmes basés sur elle. Peut-être que sa spécialisation est basée sur la matrice zéro que vous lui passez, ou quelque chose comme ça? Peut-être est-ce la mise en cache des appels de fonction répétés, qui perturbe vos timings? Peut-être qu'il optimise hors répété inutilisés produits de matrice?

pour vous prémunir contre de telles choses, utilisez une matrice de nombres aléatoires, et assurez-vous que vous forcez l'exécution en imprimant le résultat à l'écran ou sur le disque ou quelque chose comme ça.

utilisant des doubles et un tableau solide au lieu de trois pistes séparées mon code C# pour presque les mêmes résultats que C++/Java (avec votre code: 1024 - était un peu plus rapide, 2048 - était environ 140s et 4096 - était environ 22 minutes)

                1024x1024   2048x2048   4096x4096
                ---------   ---------   ---------
your C++ (ms)   6137.10     64369.29     551390.93
my C# (ms)      9730.00     90875.00    1062156.00

voici mon code:

    const int rozmer = 1024;
    double[][] matice1 = new double[rozmer * 3][];
    Random rnd = new Random();

    public Form1()
    {
        InitializeComponent();

        System.Threading.Thread thr = new System.Threading.Thread(new System.Threading.ThreadStart(() =>
        {

            string res = "";
            Stopwatch timer = new Stopwatch();
            timer.Start();

            double temp = 0;
            int r2 = rozmer * 2;

            for (int i = 0; i < rozmer*3; i++)
            {
                if (matice1[i] == null)
                {
                    matice1[i] = new double[rozmer];
                    {
                        for (int e = 0; e < rozmer; e++)
                        {
                            matice1[i][e] = rnd.NextDouble();
                        }
                    }
                }
            }
            timer.Stop();
            res += timer.ElapsedMilliseconds.ToString();

            int j = 0; int k = 0; int m = 0;

            timer.Reset();
            timer.Start();
            for (j = 0; j < rozmer; j++)
            {
                for (k = 0; k < rozmer; k++)
                {
                    temp = 0;
                    for (m = 0; m < rozmer; m++)
                    {
                        temp = temp + matice1[j][m] * matice1[m + rozmer][k];
                    }
                    matice1[j + r2][k] = temp;
                }
            }
            timer.Stop();
            this.Invoke((Action)delegate
            {
                this.Text = res + " : " + timer.ElapsedMilliseconds.ToString();
            });
        }));
        thr.Start();
    }

avez-vous vérifié que toutes les implémentations utilisaient des optimisations multi-threading pour l'algorithme ? Et ont-ils utilisé le même algorithme de multiplication ?

j'en doute.

Matlab n'est pas forcément rapide, vous avez probablement utilisé des implémentations lentes.

algorithmes pour la multiplication matricielle efficace

la réponse générale à "pourquoi matlab est plus rapide à faire xxx que les autres programmes" est que matlab a beaucoup de fonctions intégrées et optimisées.

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

Cela peut être interprété de deux façons:

1) la voie commune/théorique: Matlab n'est pas beaucoup plus rapide, vous faites juste le benchmark mal

2) la façon réaliste: pour ce truc Matlab est plus rapide dans la pratique parce que les langues comme c++ sont tout simplement trop facilement utilisées de manière inefficace.

MATLAB utilise une implémentation hautement optimisée de LAPACK à partir D'Intel connu sous le nom de Intel Math Kernel Library (Intel MKL) - en particulier la fonction dgemm . La vitesse de cette bibliothèque tire profit des fonctionnalités du processeur, y compris les instructions SIMD et les processeurs multi-core. Ils ne documentent pas quel algorithme spécifique ils utilisent. Si vous deviez appeler Intel MKL à partir de C++, vous devriez voir des performances similaires.

Je ne suis pas bien sûr ce que MATLAB bibliothèque utilise pour la multiplication GPU, mais probablement quelque chose comme NVIDIA CUBLAS .

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

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

si vous utilisez, dites, Eigen [ http://eigen.tuxfamily.org/dox/GettingStarted.html ], multithreading est intégré et le nombre de fils est réglable.

le contraste aigu n'est pas seulement dû à L'optimisation étonnante de Matlab (comme discuté par beaucoup d'autres réponses déjà), mais aussi dans la façon dont vous avez formulé matrice comme un objet.

on dirait que vous avez fait à matrix une liste de listes? Une liste de listes contient des pointeurs vers des listes qui contiennent ensuite vos éléments de matrice. Les emplacements des listes contenues sont assignés arbitrairement. Comme vous êtes en boucle sur votre premier indice (numéro de ligne?), le temps d'accès à la mémoire est très significatif. En comparaison, pourquoi ne pas essayer d'implémenter matrix en tant que 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 le titre de comparaison, vous montrera le temps d'accès à la mémoire peut être considérable!