Skip to content

Come fa numpy ad essere così veloce?

I nostri migliori programmatori hanno esaurito le loro scorte di caffè, cercando continuamente la soluzione, finché Pedro non ha trovato la soluzione in Gitea, quindi ora la condivide qui.

Soluzione:

Come sottolineato da alcuni commenti, numpy utilizza SIMD nella sua implementazione e non alloca la memoria nel punto di calcolo. Se elimino l'allocazione della memoria dalla vostra implementazione, preallocando tutti i buffer prima del calcolo, ottengo un tempo migliore rispetto a numpy anche con la versione scaler (cioè quella senza alcuna ottimizzazione).

Sempre in termini di SIMD, il motivo per cui la vostra implementazione non ha prestazioni migliori rispetto allo scaler è che i vostri schemi di accesso alla memoria non sono ideali per l'uso SIMD - fate memcopy e caricate nei registri SIMD da posizioni molto distanti l'una dall'altra - ad esempio riempite i vettori dalla riga 0 e dalla riga 511, il che potrebbe non funzionare bene con la cache o con il prefetcher SIMD.

C'è anche un errore nel modo in cui si caricano i registri SIMD (se ho capito bene quello che si sta cercando di calcolare): un registro SIMD a 256 bit può caricare 8 numeri in virgola mobile a singola precisione 8 * 32 = 256 ma nel tuo ciclo salti k per "256/sizeof(float)" che è 256/4 = 64; _x e _res sono puntatori a virgola mobile e gli intrinseci SIMD si aspettano anch'essi puntatori a virgola mobile come argomenti, quindi invece di leggere tutti gli elementi di queste righe ogni 8 virgola mobile, li leggete ogni 64 virgola mobile.

Il calcolo può essere ulteriormente ottimizzato cambiando gli schemi di accesso, ma anche osservando che si ripetono alcuni calcoli: per esempio, quando si itera con linea0 come base si calcola linea0 - linea1 ma in un momento successivo, quando si itera con linea1 come base, è necessario calcolare linea1 - linea0 che è fondamentalmente -(linea0 - linea1) cioè per ogni riga dopo riga0 molti risultati potrebbero essere riutilizzati da calcoli precedenti.
Molto spesso l'uso di SIMD o la parallelizzazione richiedono di cambiare il modo in cui si accede ai dati o si ragiona su di essi per fornire miglioramenti significativi.

Ecco quello che ho fatto come primo passo sulla base della vostra implementazione iniziale ed è più veloce di numpy (non fate caso alle cose OpenMP perché non è come dovrebbe essere fatto, volevo solo vedere come si comporta provando il modo ingenuo).

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 
#include    // compile with -mavx -msse4.1
#include 
#include 

#include 
#include 
#include 
#include 
#include 

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(stop - start);
        if (duration < best_scaler)
        {
            best_scaler = duration;
        }
    }
    cout << "Time scaler version: " << best_scaler.count() << " msn";

    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(stop - start);
     if (duration < best_simd)
    {
        best_simd = duration;
    }
}
cout << "Time SIMD version: " << best_simd.count() << " msn";

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(stop - start);
     if (duration < best_simd_2)
    {
        best_simd_2 = duration;
    }
}
cout << "Time SIMD 2 version: " << best_simd_2.count() << " msn";

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(stop - start);
     if (duration < best_simd_3)
    {
        best_simd_3 = duration;
    }
}
cout << "Time SIMD 3 version: " << best_simd_3.count() << " msn";

    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(stop - start);
         if (duration < best_par)
        {
            best_par = duration;
        }
    }
    cout << "Time OpenMP version: " << best_par.count() << " msn";

    cout << "Verificationn";
    if (equal(output, output + output_size, output_simd))
    {
        cout << "PASSEDn";
    }
    else
    {
        cout << "FAILEDn";
    }

    return 0;
}

Modifica: Piccola correzione: c'era una chiamata sbagliata relativa alla seconda versione dell'implementazione SIMD.

Come si può vedere ora, la seconda implementazione è la più veloce perché si comporta meglio dal punto di vista della località di riferimento della cache. Gli esempi 2 e 3 di implementazioni SIMD servono a illustrare come la modifica dei modelli di accesso alla memoria possa influenzare le prestazioni delle ottimizzazioni SIMD.
Per riassumere (sapendo che sono ben lungi dall'essere completo nei miei consigli), fate attenzione agli schemi di accesso alla memoria e ai carichi e ai depositi verso l'unità SIMD; la SIMD è un'unità hardware diversa all'interno del core del processore, quindi c'è una penalità nel rimescolare i dati avanti e indietro, quindi quando caricate un registro dalla memoria cercate di fare il maggior numero possibile di operazioni con quel dato e non siate troppo ansiosi di memorizzarlo di nuovo (naturalmente, nel vostro esempio potrebbe essere tutto ciò che dovete fare con i dati). Tenete inoltre presente che il numero di registri SIMD disponibili è limitato e che se ne caricate troppi, essi si "riverseranno", cioè verranno memorizzati in posizioni temporanee nella memoria principale dietro le quinte, annullando tutti i vostri guadagni. L'ottimizzazione SIMD è un vero e proprio gioco di equilibri!

C'è qualche sforzo per inserire un wrapper intrinseco multipiattaforma nello standard (ne ho sviluppato uno closed source nel mio glorioso passato) e anche se è lontano dall'essere completo, vale la pena di darci un'occhiata (leggete i documenti di accompagnamento se siete veramente interessati a imparare come funziona la SIMD).
https://github.com/VcDevel/std-simd

Questo è un complemento alla risposta postata da @celakev .
Credo di essere finalmente riuscito a capire quale fosse esattamente il problema. Il problema era non sull'allocazione della memoria nella funzione principale che esegue il calcolo.

Ciò che in realtà richiedeva tempo era accedere alla nuova memoria (fresca). Credo che il malloc restituisca pagine di memoria che sono virtuali, cioè che non corrispondono alla memoria fisica effettiva, finché non vi si accede esplicitamente. Ciò che effettivamente richiede tempo è il processo di allocazione della memoria fisica al volo (che credo sia a livello di sistema operativo) quando vi si accede nel codice della funzione.

Ecco una prova. Considerate le due seguenti funzioni banali:

#include 
#include 
#include 

float* just_alloc( size_t N ) 
{
    return (float*) aligned_alloc( 32, sizeof(float)*N );
}

void just_fill( float* _arr, size_t N ) 
{
    for (size_t i = 0; i < N; i++)
        _arr[i] = 1;
}

#define Time( code_to_benchmark, cleanup_code ) 
    do { 
        double best = 9e9; 
        for( int i = 0; i < 5; i++) { 
            struct timespec start, stop; 
            clock_gettime(CLOCK_THREAD_CPUTIME_ID, &start); 
            code_to_benchmark; 
            clock_gettime(CLOCK_THREAD_CPUTIME_ID, &stop); 
            double t = (stop.tv_sec - start.tv_sec) * 1e3 + (stop.tv_nsec - start.tv_nsec) / 1e6; 
            printf("Time[%d] = %f msn", i, t); 
            if (t < best) best = t; 
            cleanup_code; 
        } 
        printf("Best of 5 for '" #code_to_benchmark "' = %f msnn", best); 
    } while(0)

int main() 
{
    const size_t N = 512;

    Time( float* arr = just_alloc(N*N*N), free(arr) );

    float* arr = just_alloc(N*N*N);
    Time( just_fill(arr, N*N*N), ; );
    free(arr);

    return 0;
}

Ottengo i seguenti tempi, che ora illustro in dettaglio per ciascuna delle chiamate:

Time[0] = 0.000931 ms
Time[1] = 0.000540 ms
Time[2] = 0.000523 ms
Time[3] = 0.000524 ms
Time[4] = 0.000521 ms
Best of 5 for 'float* arr = just_alloc(N*N*N)' = 0.000521 ms

Time[0] = 189.822237 ms
Time[1] = 45.041083 ms
Time[2] = 46.331428 ms
Time[3] = 44.729433 ms
Time[4] = 42.241279 ms
Best of 5 for 'just_fill(arr, N*N*N)' = 42.241279 ms

Come si può vedere, l'allocazione della memoria è velocissima, ma la prima volta che si accede alla memoria è 5 volte più lenta delle altre volte. Quindi, in pratica, la ragione per cui il mio codice era lento era che ogni volta riallocavo memoria fresca che non aveva ancora un indirizzo fisico. (Correggetemi se sbaglio, ma credo che questo sia il succo del discorso).

Qui puoi vedere i commenti e le valutazioni degli utenti

Alla fine di questo articolo puoi trovare i commenti di altri creatori, puoi comunque mostrare il tuo se ti piace.



Utilizzate il nostro motore di ricerca

Ricerca
Generic filters

Lascia un commento

Il tuo indirizzo email non sarà pubblicato.