¿Por qué la multiplicación de matrices con Numpy es mucho más rápida que gsl_blas_sgemm de GSL, por ejemplo:
import numpy as np import time N = 1000 M = np.zeros(shape=(N, N), dtype=np.float) for i in range(N): for j in range(N): M[i, j] = 0.23 + 100*i + j tic = time.time() np.matmul(M, M) toc = time.time() print(toc - tic)da algo entre 0.017 - 0.019 segundos, mientras que en C++:
#include <chrono> #include <iostream> #include <gsl/gsl_matrix.h> #include <gsl/gsl_blas.h> using namespace std::chrono; int main(void) { int N = 1000; gsl_matrix_float* M = gsl_matrix_float_alloc(N, N); for (int i = 0; i < N; i++) { for (int j = 0; j < N; j++) { gsl_matrix_float_set(M, i, j, 0.23 + 100 * i + j); } } gsl_matrix_float* C = gsl_matrix_float_alloc(N, N); // save the result into C auto start = high_resolution_clock::now(); gsl_blas_sgemm(CblasNoTrans, CblasNoTrans, 1.0, M, M, 0.0, C); auto stop = high_resolution_clock::now(); auto duration = duration_cast<milliseconds>(stop - start); std::cout << duration.count() << std::endl; return 0; } Obtengo un tiempo de ejecución de la multiplicación de aproximadamente 2,7 segundos. También estoy compilando con la opción de velocidad máxima /02 . Estoy trabajando con Visual Studio. Debo hacer algo muy mal. No esperaba un rendimiento mucho mejor del código C++ porque soy consciente de que Numpy es un código C optimizado, pero tampoco esperaba que fuera unas 150 veces más lento que Python. ¿Porqué es eso? ¿Cómo puedo mejorar el tiempo de ejecución de la multiplicación en relación con Numpy?
Antecedentes del problema: necesito evaluar una integral de 1000 a 2000 dimensiones, y lo estoy haciendo con el método Monte-Carlo. Para eso, escribí casi todo el integrando como operaciones de matriz Numpy, esto funciona bastante rápido, pero lo necesito aún más rápido para evaluar el mismo integrando de 100.000 a 500.000 veces, por lo que cualquier pequeña mejora ayudaría. ¿Tiene sentido escribir el mismo código en C/C++ o debo seguir con Numpy? ¡Gracias!
TL; DR: el código C ++ y Numpy no usan la misma biblioteca de multiplicación de matrices.
La multiplicación de matrices de la biblioteca GSL no está optimizada . En mi máquina, se ejecuta secuencialmente, no usa instrucciones SIMD ( SSE / AVX ), no desenrolla eficientemente los bucles para realizar el mosaico de registros. También sospecho que tampoco usa el caché de la CPU de manera eficiente debido a la falta de mosaico. Estas optimizaciones son fundamentales para lograr un alto rendimiento y se utilizan ampliamente en bibliotecas de álgebra lineal rápida.
Numpy usa una biblioteca BLAS instalada en su máquina . En muchas plataformas Linux, utiliza OpenBLAS o Intel MKL. Ambos son muy rápidos (utilizan todos los métodos descritos anteriormente) y deberían ejecutarse en paralelo.
Puede encontrar qué implementación de BLAS utiliza Numpy aquí . En mi máquina Linux, Numpy usa por defecto CBLAS que internamente usa OpenBLAS (OpenBLAS extrañamente no es detectado directamente por Numpy).
Hay muchas implementaciones BLAS paralelas rápidas (GotoBLAS, ATLAS, BLIS, etc.). La biblioteca BLIS de código abierto es excelente porque su multiplicación de matrices es muy rápida en muchas arquitecturas diferentes.
Como resultado, la forma más sencilla de mejorar su código C++ es usar la función cblas_sgemm CBLAS y vincular una biblioteca BLAS rápida como OpenBLAS o BLIS, por ejemplo.
Para más información:
Una forma simple de ver qué tan mal funciona el GSL es usar un generador de perfiles (como perf en Linux o VTune en Windows). En su caso, rendimiento de Linux, informe que> 99% del tiempo se gasta en libgslcblas.so (es decir, la biblioteca GSL). Más específicamente, la mayor parte del tiempo de ejecución se gasta en este siguiente ciclo de ensamblaje:
250: movss (%rdx),%xmm1 add $0x4,%rax add $0x4,%rdx mulss %xmm2,%xmm1 # scalar instructions addss -0x4(%rax),%xmm1 movss %xmm1,-0x4(%rax) cmp %rax,%r9 ↑ jne 250 En cuanto a Numpy, el 99 % de su tiempo se dedica a libopenblasp-r0.3.13.so (es decir, la biblioteca OpenBLAS). Más concretamente en el siguiente código ensamblador de la función dgemm_kernel_HASWELL :
110: lea 0x80(%rsp),%rsi add $0x60,%rsi mov %r12,%rax sar $0x3,%rax cmp $0x2,%rax ↓ jl d26 prefetcht0 0x200(%rdi) # Data prefetching vmovups -0x60(%rsi),%ymm1 prefetcht0 0xa0(%rsi) vbroadcastsd -0x80(%rdi),%ymm0 # Fast SIMD instruction (AVX) prefetcht0 0xe0(%rsi) vmovups -0x40(%rsi),%ymm2 prefetcht0 0x120(%rsi) vmovups -0x20(%rsi),%ymm3 vmulpd %ymm0,%ymm1,%ymm4 prefetcht0 0x160(%rsi) vmulpd %ymm0,%ymm2,%ymm8 vmulpd %ymm0,%ymm3,%ymm12 prefetcht0 0x1a0(%rsi) vbroadcastsd -0x78(%rdi),%ymm0 vmulpd %ymm0,%ymm1,%ymm5 vmulpd %ymm0,%ymm2,%ymm9 [...]Podemos ver claramente que el código GSL no está optimizado (debido al código escalar y al bucle simple ingenuo) y que el código OpenBLAS está optimizado ya que utiliza al menos instrucciones SIMD amplias, búsqueda previa de datos y despliegue de bucle. Tenga en cuenta que el código OpenBLAS ejecutado no es óptimo, ya que podría usar las instrucciones FMA disponibles en mi procesador.