¿Por qué MATLAB es tan rápido en la multiplicación de matrices?

Estoy haciendo algunos puntos de referencia con CUDA, C ++, C #, Java y utilizando MATLAB para la verificación y la generación de matrices. Cuando realizo la multiplicación de matrices con MATLAB, 2048x2048e incluso las matrices más grandes se multiplican casi instantáneamente.

             1024x1024   2048x2048   4096x4096
             ---------   ---------   ---------
CUDA C (ms)      43.11      391.05     3407.99
C++ (ms)       6137.10    64369.29   551390.93
C# (ms)       10509.00   300684.00  2527250.00
Java (ms)      9149.90    92562.28   838357.94
MATLAB (ms)      75.01      423.10     3133.90

Solo CUDA es competitivo, pero pensé que al menos C ++ será algo cercano y no 60 veces más lento. Tampoco sé qué pensar sobre los resultados de C #. El algoritmo es el mismo que C ++ y Java, pero hay un salto gigante 2048desde 1024.

¿Cómo está MATLAB realizando la multiplicación de matrices tan rápido?

Código C ++:

float temp = 0;
timer.start();
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] * matice2[m][k];
        }
        matice3[j][k] = temp;
    }
}
timer.stop();

Aquí están mis resultados usando MATLAB R2011a + Parallel Computing Toolbox en una máquina con un Tesla C2070:

>> A = rand(1024); gA = gpuArray(A);
% warm up by executing the operations a couple of times, and then:
>> tic, C = A * A; toc
Elapsed time is 0.075396 seconds.
>> tic, gC = gA * gA; toc
Elapsed time is 0.008621 seconds.

MATLAB utiliza bibliotecas altamente optimizadas para la multiplicación de matrices, razón por la cual la multiplicación de matrices de MATLAB es tan rápida. La gpuArrayversión usa MAGMA .

Actualización mediante R2014a en una máquina con un Tesla K20c, y el nuevo timeity gputimeitfunciones:

>> A = rand(1024); gA = gpuArray(A);
>> timeit(@()A*A)
ans =
    0.0324
>> gputimeit(@()gA*gA)
ans =
    0.0022

Actualización usando R2018b en una máquina WIN64 con 16 núcleos físicos y un Tesla V100:

>> timeit(@()A*A)
ans =
    0.0229
>> gputimeit(@()gA*gA)
ans =
   4.8019e-04

(NB: en algún momento (no gpuArrayrecuerdo cuándo exactamente) cambió de MAGMA a cuBLAS, aunque MAGMA todavía se usa para algunas gpuArrayoperaciones)

Este tipo de pregunta es recurrente y debe responderse con más claridad que "MATLAB usa bibliotecas altamente optimizadas" o "MATLAB usa MKL" por una vez en Stack Overflow.

Historia:

La multiplicación de matrices (junto con la matriz de vectores, la multiplicación de vectores y muchas de las descomposiciones de matrices) es (son) los problemas más importantes en álgebra lineal. Los ingenieros han estado resolviendo estos problemas con las computadoras desde los primeros días.

No soy un experto en la historia, pero aparentemente en aquel entonces, todos reescribieron su versión FORTRAN con bucles simples. Luego llegó una cierta estandarización, con la identificación de "núcleos" (rutinas básicas) que la mayoría de los problemas de álgebra lineal necesitaban para ser resueltos. Estas operaciones básicas se estandarizaron luego en una especificación llamada: Subprogramas de álgebra lineal básica (BLAS). Los ingenieros podrían llamar a estas rutinas BLAS estándar y bien probadas en su código, haciendo su trabajo mucho más fácil.

BLAS:

BLAS evolucionó desde el nivel 1 (la primera versión que definió las operaciones escalar-vector y vector-vector) al nivel 2 (operaciones de matriz de vectores) al nivel 3 (operaciones de matriz-matriz), y proporcionó más y más "núcleos" tan estandarizados. y más de las operaciones fundamentales de álgebra lineal. Las implementaciones originales de FORTRAN 77 todavía están disponibles en el sitio web de Netlib .

Hacia un mejor rendimiento:

Entonces, a lo largo de los años (especialmente entre las versiones de nivel 1 y nivel 2 de BLAS: principios de los 80), el hardware cambió, con el advenimiento de las operaciones vectoriales y las jerarquías de caché. Estas evoluciones permitieron aumentar sustancialmente el rendimiento de las subrutinas BLAS. Luego, diferentes proveedores llegaron junto con su implementación de rutinas BLAS que eran cada vez más eficientes.

No conozco todas las implementaciones históricas (no nací o era un niño en ese entonces), pero dos de las más notables salieron a principios de la década de 2000: Intel MKL y GotoBLAS. Su Matlab utiliza el Intel MKL, que es un BLAS muy bueno y optimizado, y eso explica el gran rendimiento que ve.

Detalles técnicos sobre la multiplicación de matrices:

Entonces, ¿por qué Matlab (el MKL) es tan rápido en dgemm(multiplicación matriz-matriz general de doble precisión)? En términos simples: porque utiliza la vectorización y el buen almacenamiento en caché de datos. En términos más complejos: vea el artículo proporcionado por Jonathan Moore.

Básicamente, cuando realiza su multiplicación en el código C ++ que proporcionó, no es compatible con la caché. Dado que sospecho que creó una matriz de punteros para alinear las matrices, sus accesos en su bucle interno a la columna k-ésima de "matice2": matice2[m][k]son muy lentos. De hecho, cuando accede matice2[0][k], debe obtener el elemento k-ésimo de la matriz 0 de su matriz. Luego, en la siguiente iteración, debe acceder matice2[1][k], que es el elemento k-ésimo de otra matriz (la matriz 1). Luego, en la siguiente iteración, accede a otra matriz, y así sucesivamente ... Dado que toda la matriz matice2no puede caber en las cachés más altas (es de 8*1024*1024bytes grandes), el programa debe recuperar el elemento deseado de la memoria principal, perdiendo una gran cantidad de hora.

Si acaba de transponer la matriz, de modo que los accesos estén en direcciones de memoria contiguas, su código ya se ejecutará mucho más rápido porque ahora el compilador puede cargar filas completas en la memoria caché al mismo tiempo. Solo prueba esta versión modificada:

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();

Entonces puede ver cómo solo la localidad de caché aumentó el rendimiento de su código de manera bastante sustancial. Ahora, las dgemmimplementaciones reales explotan eso a un nivel muy extenso: realizan la multiplicación en bloques de la matriz definida por el tamaño del TLB (búfer de traducción al lado, larga historia: lo que efectivamente se puede almacenar en caché), para que se transmitan al procesador exactamente la cantidad de datos que puede procesar. El otro aspecto es la vectorización, utilizan las instrucciones vectorizadas del procesador para un rendimiento óptimo de la instrucción, lo que realmente no puede hacer desde su código C ++ multiplataforma.

Finalmente, las personas que afirman que se debe al algoritmo de Strassen o Coppersmith-Winograd están equivocadas, ambos algoritmos no son implementables en la práctica, debido a las consideraciones de hardware mencionadas anteriormente.

Esta es la razón . MATLAB no realiza una multiplicación de matriz ingenua haciendo un bucle sobre cada elemento como lo hizo en su código C ++.

Por supuesto, supongo que acabas de usar en C=A*Blugar de escribir una función de multiplicación tú mismo.

Matlab incorporó LAPACK hace algún tiempo, por lo que supongo que su multiplicación de matriz usa algo al menos tan rápido. El código fuente y la documentación de LAPACK están fácilmente disponibles.

También puede consultar el artículo de Anatomía de la multiplicación matricial de alto rendimiento de Goto y Van De Geijn en http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.140.1785&rep=rep1&type=pdf

La respuesta es LAPACK y las bibliotecas BLAS hacen que MATLAB sea cegadoramente rápido en las operaciones matriciales, no cualquier código patentado por la gente de MATLAB.

Utilice las bibliotecas LAPACK y / o BLAS en su código C ++ para operaciones matriciales y debería obtener un rendimiento similar al de MATLAB. Estas bibliotecas deberían estar disponibles gratuitamente en cualquier sistema moderno y las piezas se desarrollaron durante décadas en la academia. Tenga en cuenta que hay múltiples implementaciones, incluidas algunas de código cerrado como Intel MKL .

Una discusión sobre cómo BLAS obtiene un alto rendimiento está disponible aquí.

Por cierto, es un dolor grave en mi experiencia llamar a las bibliotecas LAPACK directamente desde c (pero vale la pena). Necesita leer la documentación MUY precisamente.

Al hacer una matriz de multiplicación, utiliza el método de multiplicación ingenuo que lleva tiempo O(n^3).

Existe un algoritmo de multiplicación de matrices que toma O(n^2.4). Lo que significa que en n=2000su algoritmo requiere ~ 100 veces más cómputo que el mejor algoritmo.
Realmente debería consultar la página de Wikipedia para la multiplicación de matrices para obtener más información sobre las formas eficientes de implementarlo.

Dependiendo de su versión de Matlab, creo que ya podría estar usando su GPU.

Otra cosa; Matlab realiza un seguimiento de muchas propiedades de su matriz; ya sea diagonal, hermetiana, etc., y especializa sus algoritmos basados ​​en ella. ¿Tal vez se especializa en base a la matriz cero que está pasando, o algo así? Tal vez es el almacenamiento en caché de llamadas de función repetidas, lo que arruina sus tiempos? ¿Quizás optimiza los productos de matriz repetidos no utilizados?

Para evitar que sucedan estas cosas, use una matriz de números aleatorios y asegúrese de forzar la ejecución imprimiendo el resultado en la pantalla, el disco o algo así.

MATLAB utiliza una implementación altamente optimizada de LAPACK de Intel conocida como Intel Math Kernel Library (Intel MKL), específicamente la función dgemm . La velocidad Esta biblioteca aprovecha las características del procesador, incluidas las instrucciones SIMD y los procesadores multinúcleo. No documentan qué algoritmo específico usan. Si llamara a Intel MKL desde C ++, debería ver un rendimiento similar.

No estoy seguro de qué biblioteca utiliza MATLAB para la multiplicación de GPU, pero probablemente algo así como nVidia CUBLAS .

La respuesta general a "¿Por qué Matlab es más rápido para hacer xxx que otros programas?" Es que Matlab tiene muchas funciones integradas y optimizadas.

Los otros programas que se usan a menudo no tienen estas funciones, por lo que las personas aplican sus propias soluciones creativas, que son sorprendentemente más lentas que el código optimizado profesionalmente.

Esto se puede interpretar de dos maneras:

1) La forma común / teórica: Matlab no es significativamente más rápido, solo está haciendo mal el punto de referencia

2) La forma realista: para estas cosas, Matlab es más rápido en la práctica porque los lenguajes como c ++ se usan con demasiada facilidad de manera ineficaz.

El fuerte contraste no solo se debe a la asombrosa optimización de Matlab (como ya se discutió en muchas otras respuestas), sino también a la forma en que formuló la matriz como un objeto.

Parece que hiciste matriz una lista de listas? Una lista de listas contiene punteros a listas que luego contienen sus elementos de matriz. Las ubicaciones de las listas contenidas se asignan arbitrariamente. A medida que recorre su primer índice (¿número de fila?), El tiempo de acceso a la memoria es muy significativo. En comparación, ¿por qué no intentas implementar la matriz como una sola lista / vector usando el siguiente método?

#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);
};

Y

double &matrix::operator()(int i, int j) {
    return M[n_col * i + j];
}

Se debe utilizar el mismo algoritmo de multiplicación para que el número de flop sea el mismo. (n ^ 3 para matrices cuadradas de tamaño n)

Le pido que cronometre para que el resultado sea comparable al que tenía anteriormente (en la misma máquina). ¡Con la comparación, mostrará exactamente cuán significativo puede ser el tiempo de acceso a la memoria!

Es lento en C ++ porque no estás usando subprocesos múltiples. Esencialmente, si A = BC, donde todas son matrices, la primera fila de A puede calcularse independientemente de la segunda fila, etc. Si A, B y C son todas matrices n por n, puede acelerar la multiplicación por un factor de n ^ 2, como

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

Si usa, por ejemplo, Eigen [ http://eigen.tuxfamily.org/dox/GettingStarted.html ], el subprocesamiento múltiple está integrado y el número de subprocesos es ajustable.

Porque MATLAB es un lenguaje de programación desarrollado inicialmente para álgebra lineal numérica (manipulaciones matriciales), que tiene bibliotecas especialmente desarrolladas para multiplicaciones matriciales. Y ahora MATLAB también puede usar las GPU (unidad de procesamiento de gráficos) para esto adicionalmente.

Y si miramos los resultados de su cálculo:

             1024x1024   2048x2048   4096x4096
             ---------   ---------   ---------
CUDA C (ms)      43.11      391.05     3407.99
C++ (ms)       6137.10    64369.29   551390.93
C# (ms)       10509.00   300684.00  2527250.00
Java (ms)      9149.90    92562.28   838357.94
MATLAB (ms)      75.01      423.10     3133.90

entonces podemos ver que no solo MATLAB es tan rápido en la multiplicación de matrices: CUDA C (lenguaje de programación de NVIDIA) tiene mejores resultados que MATLAB. CUDA C también tiene bibliotecas especialmente desarrolladas para multiplicaciones matriciales y utiliza las GPU.

Breve historia de MATLAB

Cleve Moler, presidente del departamento de informática de la Universidad de Nuevo México, comenzó a desarrollar MATLAB a fines de la década de 1970. Lo diseñó para dar a sus alumnos acceso a LINPACK (una biblioteca de software para realizar álgebra lineal numérica) y EISPACK(es una biblioteca de software para el cálculo numérico de álgebra lineal) sin que tengan que aprender Fortran. Pronto se extendió a otras universidades y encontró una audiencia fuerte dentro de la comunidad de matemáticas aplicadas. Jack Little, un ingeniero, estuvo expuesto durante una visita que Moler realizó a la Universidad de Stanford en 1983. Reconociendo su potencial comercial, se unió a Moler y Steve Bangert. Reescribieron MATLAB en C y fundaron MathWorks en 1984 para continuar su desarrollo. Estas bibliotecas reescritas se conocían como JACKPAC. En 2000, MATLAB fue reescrito para usar un nuevo conjunto de bibliotecas para la manipulación de matrices, LAPACK (es una biblioteca de software estándar para álgebra lineal numérica).

Fuente

¿Qué es CUDA C?

CUDA C también utiliza bibliotecas especialmente desarrolladas para multiplicaciones matriciales como OpenGL (Open Graphics Library). Utiliza también GPU y Direct3D (en MS Windows).

La plataforma CUDA está diseñada para trabajar con lenguajes de programación como C, C ++ y Fortran. Esta accesibilidad facilita a los especialistas en programación paralela el uso de recursos de GPU, en contraste con API anteriores como Direct3D y OpenGL , que requerían habilidades avanzadas en programación de gráficos. Además, CUDA admite marcos de programación como OpenACC y OpenCL .

Ejemplo de flujo de procesamiento de CUDA:

1. Copie datos de la memoria principal a la memoria de la GPU
2. CPU inicia el núcleo de cómputo GPU
3. Los núcleos CUDA de GPU ejecutan el núcleo en paralelo
4. Copie los datos resultantes de la memoria de la GPU a la memoria principal

Comparación de velocidades de ejecución de CPU y GPU

Ejecutamos un punto de referencia en el que medimos la cantidad de tiempo que llevó ejecutar 50 pasos de tiempo para tamaños de cuadrícula de 64, 128, 512, 1024 y 2048 en un procesador Intel Xeon X5650 y luego usando una GPU NVIDIA Tesla C2050.

Para un tamaño de cuadrícula de 2048, el algoritmo muestra una disminución de 7,5 veces en el tiempo de cálculo de más de un minuto en la CPU a menos de 10 segundos en la GPU. El diagrama de escala de registro muestra que la CPU es realmente más rápida para tamaños de cuadrícula pequeños. Sin embargo, a medida que la tecnología evoluciona y madura, las soluciones de GPU son cada vez más capaces de manejar problemas más pequeños, una tendencia que esperamos continuar.

Fuente

De la introducción de la Guía de programación de CUDA C:

Impulsado por la insaciable demanda del mercado de gráficos 3D de alta definición en tiempo real, la Unidad de Procesador Gráfico o GPU programable ha evolucionado hasta convertirse en un procesador de muchos núcleos altamente paralelo, multiproceso, con una enorme potencia computacional y un ancho de banda de memoria muy alto, como se ilustra por Figure 1y Figure 2.

Figura 1. Operaciones de punto flotante por segundo para la CPU y la GPU

La Figura 2 . Ancho de banda de memoria para la CPU y GPU

La razón detrás de la discrepancia en la capacidad de punto flotante entre la CPU y la GPU es que la GPU está especializada para el cómputo altamente paralelo, intensivo en cómputo, exactamente de lo que se trata el renderizado de gráficos, y por lo tanto diseñado de tal manera que más transistores se dediquen al procesamiento de datos en lugar de almacenamiento en caché de datos y control de flujo, como se ilustra esquemáticamente por Figure 3.

Figura 3 . La GPU dedica más transistores al procesamiento de datos

Más específicamente, la GPU es especialmente adecuada para abordar problemas que pueden expresarse como cálculos de datos paralelos: el mismo programa se ejecuta en muchos elementos de datos en paralelo, con alta intensidad aritmética, la relación entre operaciones aritméticas y operaciones de memoria. Debido a que el mismo programa se ejecuta para cada elemento de datos, existe un requisito menor para un control de flujo sofisticado, y debido a que se ejecuta en muchos elementos de datos y tiene una alta intensidad aritmética, la latencia de acceso a la memoria se puede ocultar con cálculos en lugar de grandes cachés de datos .

El procesamiento paralelo de datos asigna elementos de datos a subprocesos de procesamiento paralelo. Muchas aplicaciones que procesan grandes conjuntos de datos pueden usar un modelo de programación de datos paralelos para acelerar los cálculos. En la representación 3D, grandes conjuntos de píxeles y vértices se asignan a hilos paralelos. Del mismo modo, las aplicaciones de procesamiento de imágenes y medios, como el procesamiento posterior de imágenes renderizadas, la codificación y decodificación de video, el escalado de imágenes, la visión estéreo y el reconocimiento de patrones pueden asignar bloques de imagen y píxeles a hilos de procesamiento paralelos. De hecho, muchos algoritmos fuera del campo de procesamiento y procesamiento de imágenes se aceleran mediante el procesamiento paralelo de datos, desde el procesamiento general de señales o la simulación física hasta las finanzas computacionales o la biología computacional.

Fuente

Lectura avanzada

• GPU (unidad de procesamiento de gráficos)
• MATLAB
• Guía de programación de CUDA C
• Uso de GPU en MATLAB

• Subprogramas de álgebra lineal básica (BLAS)

• Anatomía de la multiplicación matricial de alto rendimiento , de Kazushige Goto y Robert A. Van De Geijn

Algunos facs interesantes

He escrito la multiplicación de matrices de C ++ que es tan rápida como la de Matlab pero tuvo mucho cuidado. (Antes de que Matlab usara GPU para esto).

Cita de esta respuesta .