Comment Numpy est-il si rapide?

Comme le souligne certains commentaires, Numpy utilise SIMD dans son implémentation et il n'alloue pas de la mémoire au point de calcul. Si j'élimine l'allocation de mémoire de votre implémentation, pré-allocation de tous les tampons avant le calcul, je reçois un meilleur temps par rapport à Numpy même avec la version Scaller (c'est-à-dire celle sans optimisations).

également en termes de SIMD et pourquoi votre implémentation ne fonctionne pas beaucoup mieux que l'échelle est parce que vos modèles d'accès à la mémoire ne sont pas idéaux pour l'utilisation SIMD - vous faites de la mémoire et que vous chargez dans des registres SIMD à partir d'emplacements qui sont loin les uns des autres - par exemple Vous remplissez les vecteurs de la ligne 0 et de la ligne 511, qui pourraient ne pas bien jouer avec le cache ou avec le préfetcher SIMD.

Il y a aussi une erreur dans la façon dont vous chargez les registres SIMD (si je comprenais correctement ce que vous 're essaie de calculer): un registre SIMD 256 bits peut charger 8 numéros de point flottante à un seul précis (float) " qui est 256/4 = 64 ; _x et _res sont des pointeurs flottants et les intrinsèques SIMD s'attendent également à des pointeurs flottants comme arguments, donc au lieu de lire tous les éléments de ces lignes tous les 8 flotteurs que vous les lisez tous les 64 flotteurs.

Le calcul peut être optimisé davantage en modifiant les modèles d'accès mais aussi en observant que vous répétez certains calculs: par exemple Lorsque vous en itérant avec line0 en tant que base, vous calculez line0 - line1 mais à un moment futur, lorsque vous en itérant avec line1 comme base, vous devez calcul line1 - line0 qui est essentiellement - (line0 - line1) , c'est-à-dire pour chaque ligne après line0 beaucoup de résultats pourraient être réutilisés à partir de calculs précédents. Souvent, l'utilisation ou la parallélisation SIMD nécessite de modifier la façon dont les données sont accessibles ou motivées afin de fournir des améliorations significatives.

Voici ce que j'ai fait en première étape en fonction de votre implémentation initiale et informatique est plus rapide que le Numpy (ne vous dérange pas les trucs OpenMP car ce n'est pas comme ça que c'est censé être fait, je voulais juste voir comment il se comporte en essayant la manière naïve).

C++
Time scaler version: 55 ms
Time SIMD version: 53 ms
**Time SIMD 2 version: 33 ms**
Time SIMD 3 version: 168 ms
Time OpenMP version: 59 ms

Python numpy
>> best of 5 = 88.794 ms


#include <cstdlib>
#include <xmmintrin.h>   // compile with -mavx -msse4.1
#include <pmmintrin.h>
#include <immintrin.h>

#include <numeric>
#include <algorithm>
#include <chrono>
#include <iostream>
#include <cstring>

using namespace std;

float* pairwise_sub_naive (const float* input, float* output, int n) 
{
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            for (int k = 0; k < n; k++)
                output[(i * n + j) * n + k] = input[i * n + k] - input[j * n + k];
          }
    }
    return output;
}

float* pairwise_sub_simd (const float* input, float* output, int n) 
{    
    for (int i = 0; i < n; i++) 
    {
        const int idxi = i * n;
        for (int j = 0; j < n; j++)
        {
            const int idxj = j * n;
            const int outidx = idxi + j;
            for (int k = 0; k < n; k += 8) 
            {
                __m256 A = _mm256_load_ps(input + idxi + k);
                __m256 B = _mm256_load_ps(input + idxj + k);
                _mm256_store_ps(output + outidx * n + k, _mm256_sub_ps( A, B ));
            }
        }
    }
    
    return output;
}

float* pairwise_sub_simd_2 (const float* input, float* output, int n) 
{
    float* line_buffer = (float*) aligned_alloc(32, n * sizeof(float));

    for (int i = 0; i < n; i++) 
    {
        const int idxi = i * n;
        for (int j = 0; j < n; j++)
        {
            const int idxj = j * n;
            const int outidx = idxi + j;
            for (int k = 0; k < n; k += 8) 
            {
                __m256 A = _mm256_load_ps(input + idxi + k);
                __m256 B = _mm256_load_ps(input + idxj + k);
                _mm256_store_ps(line_buffer + k, _mm256_sub_ps( A, B ));
            }
            memcpy(output + outidx * n, line_buffer, n);
        }
    }
    
    return output;
}

float* pairwise_sub_simd_3 (const float* input, float* output, int n) 
{    
    for (int i = 0; i < n; i++) 
    {
        const int idxi = i * n;
        for (int k = 0; k < n; k += 8) 
        {
            __m256 A = _mm256_load_ps(input + idxi + k);
            for (int j = 0; j < n; j++)
            {
                const int idxj = j * n;
                const int outidx = (idxi + j) * n;
                __m256 B = _mm256_load_ps(input + idxj + k);
                _mm256_store_ps(output + outidx + k, _mm256_sub_ps( A, B     ));
             }
        }
    }

    return output;
}

float* pairwise_sub_openmp (const float* input, float* output, int n)
{
    int i, j;
    #pragma omp parallel for private(j)
    for (i = 0; i < n; i++) 
    {
        for (j = 0; j < n; j++)
        {
            const int idxi = i * n; 
            const int idxj = j * n;
            const int outidx = idxi + j;
            for (int k = 0; k < n; k += 8) 
            {
                __m256 A = _mm256_load_ps(input + idxi + k);
                __m256 B = _mm256_load_ps(input + idxj + k);
                _mm256_store_ps(output + outidx * n + k, _mm256_sub_ps( A, B ));
            }
        }
    }
    /*for (i = 0; i < n; i++) 
    {
        for (j = 0; j < n; j++) 
        {
            for (int k = 0; k < n; k++)
            {
                output[(i * n + j) * n + k] = input[i * n + k] - input[j * n + k];
            }
        }
    }*/
    
    return output;
}

int main ()
{
    constexpr size_t n = 512;
    constexpr size_t input_size = n * n;
    constexpr size_t output_size = n * n * n;

    float* input = (float*) aligned_alloc(32, input_size * sizeof(float));
    float* output = (float*) aligned_alloc(32, output_size * sizeof(float));

    float* input_simd = (float*) aligned_alloc(32, input_size * sizeof(float));
    float* output_simd = (float*) aligned_alloc(32, output_size * sizeof(float));

    float* input_par = (float*) aligned_alloc(32, input_size * sizeof(float));
    float* output_par = (float*) aligned_alloc(32, output_size * sizeof(float));

    iota(input, input + input_size, float(0.0));
    fill(output, output + output_size, float(0.0));

    iota(input_simd, input_simd + input_size, float(0.0));
    fill(output_simd, output_simd + output_size, float(0.0));
    
    iota(input_par, input_par + input_size, float(0.0));
    fill(output_par, output_par + output_size, float(0.0));

    std::chrono::milliseconds best_scaler{100000};
    for (int i = 0; i < 5; ++i)
    {
        auto start = chrono::high_resolution_clock::now();
        pairwise_sub_naive(input, output, n);
        auto stop = chrono::high_resolution_clock::now();

        auto duration = chrono::duration_cast<chrono::milliseconds>(stop - start);
        if (duration < best_scaler)
        {
            best_scaler = duration;
        }
    }
    cout << "Time scaler version: " << best_scaler.count() << " ms\n";

    std::chrono::milliseconds best_simd{100000};
for (int i = 0; i < 5; ++i)
{
    auto start = chrono::high_resolution_clock::now();
    pairwise_sub_simd(input_simd, output_simd, n);
    auto stop = chrono::high_resolution_clock::now();

    auto duration = chrono::duration_cast<chrono::milliseconds>(stop - start);
     if (duration < best_simd)
    {
        best_simd = duration;
    }
}
cout << "Time SIMD version: " << best_simd.count() << " ms\n";

std::chrono::milliseconds best_simd_2{100000};
for (int i = 0; i < 5; ++i)
{
    auto start = chrono::high_resolution_clock::now();
    pairwise_sub_simd_2(input_simd, output_simd, n);
    auto stop = chrono::high_resolution_clock::now();

    auto duration = chrono::duration_cast<chrono::milliseconds>(stop - start);
     if (duration < best_simd_2)
    {
        best_simd_2 = duration;
    }
}
cout << "Time SIMD 2 version: " << best_simd_2.count() << " ms\n";

std::chrono::milliseconds best_simd_3{100000};
for (int i = 0; i < 5; ++i)
{
    auto start = chrono::high_resolution_clock::now();
    pairwise_sub_simd_3(input_simd, output_simd, n);
    auto stop = chrono::high_resolution_clock::now();

    auto duration = chrono::duration_cast<chrono::milliseconds>(stop - start);
     if (duration < best_simd_3)
    {
        best_simd_3 = duration;
    }
}
cout << "Time SIMD 3 version: " << best_simd_3.count() << " ms\n";

    std::chrono::milliseconds best_par{100000};
    for (int i = 0; i < 5; ++i)
    {
        auto start = chrono::high_resolution_clock::now();
        pairwise_sub_openmp(input_par, output_par, n);
        auto stop = chrono::high_resolution_clock::now();

        auto duration = chrono::duration_cast<chrono::milliseconds>(stop - start);
         if (duration < best_par)
        {
            best_par = duration;
        }
    }
    cout << "Time OpenMP version: " << best_par.count() << " ms\n";

    cout << "Verification\n";
    if (equal(output, output + output_size, output_simd))
    {
        cout << "PASSED\n";
    }
    else
    {
        cout << "FAILED\n";
    }

    return 0;
}

petite correction car il y avait un mauvais appel lié à la deuxième version de l'implémentation SIMD.

Comme vous pouvez le voir maintenant, la deuxième implémentation est la plus rapide telle qu'elle se comporte Le meilleur du point de vue de la localité de référence du cache. Les exemples 2 et 3 des implémentations SIMD sont là pour illustrer pour vous comment changer les modèles d'accès à la mémoire pour influencer les performances de vos optimisations SIMD. Pour résumer (sachant que je suis loin d'être complet à mon conseil), soyez conscient de vos modèles d'accès à la mémoire et des charges et des magasins à \ de l'unité SIMD; Le SIMD est une unité matérielle différente à l'intérieur du noyau du processeur, il y a donc une pénalité dans la mélange de données dans les deux sens, donc lorsque vous chargez un registre à partir de la mémoire, essayez de faire autant d'opérations que possible avec ces données et ne soyez pas trop désireux de stocker Il est de retour (bien sûr, dans votre exemple, cela pourrait être tout ce que vous devez faire avec les données). Soyez conscient également qu'il y a un nombre limité de registres SIMD disponibles et si vous en chargez trop, ils "déversent", c'est-à-dire qu'ils seront retrouvés dans des emplacements temporaires en mémoire principale dans les coulisses en tuant tous vos gains. Optimisation SIMD, c'est un véritable acte d'équilibre!

Il y a des efforts pour mettre un emballage intrinsèque multiplateforme dans la norme (je me suis développé une source fermée dans mon passé glorieux) et même c'est loin d'être Complet, cela vaut la peine de jeter un œil (lisez les articles d'accompagnement si vous êtes vraiment intéressé à savoir comment fonctionne SIMD). https://github.com/vcdevel/std-simd

p >

Un peu en retard à la partie, mais je voulais ajouter une méthode paire avec eigen , qui est censé donner à C ++ une capacité de manipulation d'algèbre de haut niveau et utiliser SIMD sous le capot. Tout comme Numpy.

Voici l'implémentation

// inside the pairwise method
        for (int k = 0; k < input.cols(); ++k)
                output[k] -= input.col(k) * Eigen::RowVectorXf::Ones(input.cols());


// at allocation time
        std::vector<Eigen::MatrixXf> output(n, input);

Les tests de référence complets suggérés par @celavek sur mon système sont

Time scaler version: 57 ms
Time SIMD version: 58 ms
Time SIMD 2 version: 40 ms
Time SIMD 3 version: 58 ms
Time OpenMP version: 58 ms

Time Eigen version: 76 ms

Numpy >> best of 5 = 118.489 ms