Estoy tratando de entender cómo numpy puede ser tan rápido, según mi impactante comparación con el código C/C++ optimizado que aún está lejos de reproducir la velocidad de numpy.
Considere el siguiente ejemplo: dada una matriz 2D con shape=(N, N) y dtype=float32 , que representa una lista de N vectores de N dimensiones, estoy calculando las diferencias por pares entre cada par de vectores. Usando la transmisión numpy , esto simplemente se escribe como:
def pairwise_sub_numpy( X ): return X - X[:, None, :] Usando timeit , puedo medir el rendimiento para N=512 : toma 88 ms por llamada en mi computadora portátil.
Ahora, en C/C++ una implementación ingenua se escribe como:
#define X(i, j) _X[(i)*N + (j)] #define res(i, j, k) _res[((i)*N + (j))*N + (k)] float* pairwise_sub_naive( const float* _X, int N ) { float* _res = (float*) aligned_alloc( 32, N*N*N*sizeof(float)); for (int i = 0; i < N; i++) { for (int j = 0; j < N; j++) { for (int k = 0; k < N; k++) res(i,j,k) = X(i,k) - X(j,k); } } return _res; } Compilando usando gcc -O3 con el indicador -O3, obtengo 195 ms por llamada para pairwise_sub_naive(X) , que no está tan mal dada la simplicidad del código, pero aproximadamente 2 veces más lento que numpy .
Ahora empiezo a ponerme serio y agrego algunas pequeñas optimizaciones, indexando los vectores de fila directamente:
float* pairwise_sub_better( const float* _X, int N ) { float* _res = (float*) aligned_alloc( 32, N*N*N*sizeof(float)); for (int i = 0; i < N; i++) { const float* xi = & X(i,0); for (int j = 0; j < N; j++) { const float* xj = & X(j,0); float* r = &res(i,j,0); for (int k = 0; k < N; k++) r[k] = xi[k] - xj[k]; } } return _res; }La velocidad se mantiene igual en 195 ms, lo que significa que el compilador pudo calcular esa cantidad. Ahora usemos instrucciones vectoriales SIMD:
float* pairwise_sub_simd( const float* _X, int N ) { float* _res = (float*) aligned_alloc( 32, N*N*N*sizeof(float)); // create caches for row vectors which are memory-aligned float* xi = (float*)aligned_alloc(32, N * sizeof(float)); float* xj = (float*)aligned_alloc(32, N * sizeof(float)); for (int i = 0; i < N; i++) { memcpy(xi, & X(i,0), N*sizeof(float)); for (int j = 0; j < N; j++) { memcpy(xj, & X(j,0), N*sizeof(float)); float* r = &res(i,j,0); for (int k = 0; k < N; k += 256/sizeof(float)) { const __m256 A = _mm256_load_ps(xi+k); const __m256 B = _mm256_load_ps(xj+k); _mm256_store_ps(r+k, _mm256_sub_ps( A, B )); } } } free(xi); free(xj); return _res; }Esto solo produce un pequeño impulso (178 ms en lugar de 194 ms por llamada de función).
Entonces me preguntaba si un enfoque "en bloques", como el que se usa para optimizar los productos punto, podría ser beneficioso:
float* pairwise_sub_blocks( const float* _X, int N ) { float* _res = (float*) aligned_alloc( 32, N*N*N*sizeof(float)); #define B 8 float cache1[B*B], cache2[B*B]; for (int bi = 0; bi < N; bi+=B) for (int bj = 0; bj < N; bj+=B) for (int bk = 0; bk < N; bk+=B) { // load first 8x8 block in the cache for (int i = 0; i < B; i++) for (int k = 0; k < B; k++) cache1[B*i + k] = X(bi+i, bk+k); // load second 8x8 block in the cache for (int j = 0; j < B; j++) for (int k = 0; k < B; k++) cache2[B*j + k] = X(bj+j, bk+k); // compute local operations on the caches for (int i = 0; i < B; i++) for (int j = 0; j < B; j++) for (int k = 0; k < B; k++) res(bi+i,bj+j,bk+k) = cache1[B*i + k] - cache2[B*j + k]; } return _res; }Y, sorprendentemente, este es el método más lento hasta el momento (258 ms por llamada de función).
Para resumir, a pesar de algunos esfuerzos con un código C++ optimizado, no puedo acercarme a los 88 ms/llamada que numpy logra sin esfuerzo. ¿Alguna idea de por qué?
Nota: Por cierto, estoy deshabilitando los subprocesos múltiples numpy y, de todos modos, este tipo de operación no es multiproceso.
Editar : código exacto para comparar el código numpy:
import numpy as np def pairwise_sub_numpy( X ): return X - X[:, None, :] N = 512 X = np.random.rand(N,N).astype(np.float32) import timeit times = timeit.repeat('pairwise_sub_numpy( X )', globals=globals(), number=1, repeat=5) print(f">> best of 5 = {1000*min(times):.3f} ms")Punto de referencia completo para el código C:
#include <stdio.h> #include <string.h> #include <xmmintrin.h> // compile with -mavx -msse4.1 #include <pmmintrin.h> #include <immintrin.h> #include <time.h> #define X(i, j) _x[(i)*N + (j)] #define res(i, j, k) _res[((i)*N + (j))*N + (k)] float* pairwise_sub_naive( const float* _x, int N ) { float* _res = (float*) aligned_alloc( 32, N*N*N*sizeof(float)); for (int i = 0; i < N; i++) { for (int j = 0; j < N; j++) { for (int k = 0; k < N; k++) res(i,j,k) = X(i,k) - X(j,k); } } return _res; } float* pairwise_sub_better( const float* _x, int N ) { float* _res = (float*) aligned_alloc( 32, N*N*N*sizeof(float)); for (int i = 0; i < N; i++) { const float* xi = & X(i,0); for (int j = 0; j < N; j++) { const float* xj = & X(j,0); float* r = &res(i,j,0); for (int k = 0; k < N; k++) r[k] = xi[k] - xj[k]; } } return _res; } float* pairwise_sub_simd( const float* _x, int N ) { float* _res = (float*) aligned_alloc( 32, N*N*N*sizeof(float)); // create caches for row vectors which are memory-aligned float* xi = (float*)aligned_alloc(32, N * sizeof(float)); float* xj = (float*)aligned_alloc(32, N * sizeof(float)); for (int i = 0; i < N; i++) { memcpy(xi, & X(i,0), N*sizeof(float)); for (int j = 0; j < N; j++) { memcpy(xj, & X(j,0), N*sizeof(float)); float* r = &res(i,j,0); for (int k = 0; k < N; k += 256/sizeof(float)) { const __m256 A = _mm256_load_ps(xi+k); const __m256 B = _mm256_load_ps(xj+k); _mm256_store_ps(r+k, _mm256_sub_ps( A, B )); } } } free(xi); free(xj); return _res; } float* pairwise_sub_blocks( const float* _x, int N ) { float* _res = (float*) aligned_alloc( 32, N*N*N*sizeof(float)); #define B 8 float cache1[B*B], cache2[B*B]; for (int bi = 0; bi < N; bi+=B) for (int bj = 0; bj < N; bj+=B) for (int bk = 0; bk < N; bk+=B) { // load first 8x8 block in the cache for (int i = 0; i < B; i++) for (int k = 0; k < B; k++) cache1[B*i + k] = X(bi+i, bk+k); // load second 8x8 block in the cache for (int j = 0; j < B; j++) for (int k = 0; k < B; k++) cache2[B*j + k] = X(bj+j, bk+k); // compute local operations on the caches for (int i = 0; i < B; i++) for (int j = 0; j < B; j++) for (int k = 0; k < B; k++) res(bi+i,bj+j,bk+k) = cache1[B*i + k] - cache2[B*j + k]; } return _res; } int main() { const int N = 512; float* _x = (float*) malloc( N * N * sizeof(float) ); for( int i = 0; i < N; i++) for( int j = 0; j < N; j++) X(i,j) = ((i+j*j+17*i+101) % N) / float(N); double best = 9e9; for( int i = 0; i < 5; i++) { struct timespec start, stop; clock_gettime(CLOCK_THREAD_CPUTIME_ID, &start); //float* res = pairwise_sub_naive( _x, N ); //float* res = pairwise_sub_better( _x, N ); //float* res = pairwise_sub_simd( _x, N ); float* res = pairwise_sub_blocks( _x, N ); clock_gettime(CLOCK_THREAD_CPUTIME_ID, &stop); double t = (stop.tv_sec - start.tv_sec) * 1e6 + (stop.tv_nsec - start.tv_nsec) / 1e3; // in microseconds if (t < best) best = t; free( res ); } printf("Best of 5 = %f ms\n", best / 1000); free( _x ); return 0; } Compilado usando gcc 7.3.0 gcc -Wall -O3 -mavx -msse4.1 -o test_simd test_simd.c
Resumen de tiempos en mi máquina:
| Implementación | Hora |
|---|---|
| entumecido | 88ms |
| C++ ingenuo | 194ms |
| C++ mejor | 195ms |
| C++ SIMD | 178ms |
| C++ bloqueado | 258ms |
| C++ bloqueado (gcc 8.3.1) | 217ms |
Este es un complemento a la respuesta publicada por @celakev. Creo que finalmente llegué a entender cuál era exactamente el problema. El problema no se trataba de asignar la memoria en la función principal que realiza el cálculo.
Lo que realmente tomaba tiempo era acceder a la memoria nueva (nueva) . Creo que la llamada malloc devuelve páginas de memoria que son virtuales, es decir, que no corresponden a la memoria física real, hasta que se accede explícitamente. Lo que realmente lleva tiempo es el proceso de asignación de memoria física sobre la marcha (que creo que es a nivel del sistema operativo) cuando se accede a ella en el código de función.
Aquí hay una prueba. Considere las dos funciones triviales siguientes:
#include <stdio.h> #include <stdlib.h> #include <time.h> 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 ms\n", i, t); \ if (t < best) best = t; \ cleanup_code; \ } \ printf("Best of 5 for '" #code_to_benchmark "' = %f ms\n\n", 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; }Obtengo los siguientes tiempos, que ahora detallo para cada una de las llamadas:
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 msComo puede ver, la asignación de memoria es increíblemente rápida, pero la primera vez que se accede a la memoria, es 5 veces más lenta que las otras veces. Entonces, básicamente, la razón por la que mi código era lento era porque cada vez estaba reasignando memoria nueva que aún no tenía una dirección física. (Corrígeme si me equivoco, ¡pero creo que esa es la esencia!)
Un poco tarde para la fiesta, pero quería agregar un método por pairwise con Eigen , que se supone que le da a C++ una capacidad de manipulación de álgebra de alto nivel y usa SIMD bajo el capó. Al igual que numpy.
Aquí está la implementación
#include <iostream> #include <vector> #include <chrono> #include <algorithm> #include <Eigen/Dense> auto pairwise_eigen(const Eigen::MatrixXf &input, std::vector<Eigen::MatrixXf> &output) { for (int k = 0; k < input.cols(); ++k) output[k] = input // subtract matrix with repeated k-th column - input.col(k) * Eigen::RowVectorXf::Ones(input.cols()); } int main() { constexpr size_t n = 512; // allocate input and output Eigen::MatrixXf input = Eigen::MatrixXf::Random(n, n); std::vector<Eigen::MatrixXf> output(n); std::chrono::milliseconds best_eigen{100000}; for (int i = 0; i < 5; ++i) { auto start = std::chrono::high_resolution_clock::now(); pairwise_eigen(input, output); auto end = std::chrono::high_resolution_clock::now(); auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(end-start); if (duration < best_eigen) best_eigen = duration; } std::cout << "Time Eigen version: " << best_eigen.count() << " ms\n"; return 0; }Las pruebas de referencia completas sugeridas por @celavek en mi sistema son
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 msCon Eigen todavía hay una mejora notable con respecto a Numpy, pero no tan impresionante en comparación con las implementaciones "en bruto" (ciertamente hay algo de sobrecarga). Una optimización adicional es asignar el vector de salida con copias de la entrada y luego restar directamente de cada entrada de vector, simplemente reemplazando las siguientes líneas
// 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);Esto empuja lo mejor de 5 a 60 ms.