ejemplo - function matlab

¿por qué a*b*a toma más tiempo que(a ''*(a*b)'') ''cuando se usa gpuArray en scripts de Matlab?

Question

El siguiente código realiza la operación de la misma operación en gpuArrays a y b de dos maneras diferentes. La primera parte calcula (a''*(a*b)'')'' , mientras que la segunda parte calcula a*b*a . Luego se verifica que los resultados son los mismos.

%function test
clear
rng(''default'');rng(1);
a=sprand(3000,3000,0.1);
b=rand(3000,3000);
a=gpuArray(a);
b=gpuArray(b);
tic;
c1=gather(transpose(transpose(a)*transpose(a*b)));
disp([''time for (a''''*(a*b)'''')'''': '' , num2str(toc),''s''])

clearvars -except c1

rng(''default'');
rng(1)
a=sprand(3000,3000,0.1);
b=rand(3000,3000);
a=gpuArray(a);
b=gpuArray(b);
tic;
c2=gather(a*b*a);
disp([''time for a*b*a: '' , num2str(toc),''s''])

disp([''error = '',num2str(max(max(abs(c1-c2))))])

%end

Sin embargo, computar (a''*(a*b)'')'' es aproximadamente 4 veces más rápido que computar a*b*a . Aquí está la salida de la secuencia de comandos anterior en R2018a en un Nvidia K20 (he probado diferentes versiones y diferentes GPU con el mismo comportamiento).

>> test
time for (a''*(a*b)'')'': 0.43234s
time for a*b*a: 1.7175s
error = 2.0009e-11

Aún más extraño, si las primeras y últimas líneas del script anterior no tienen comentarios (para convertirlo en una función), entonces ambas toman la mayor cantidad de tiempo (~ 1.7s en lugar de ~ 0.4s). A continuación se muestra la salida para este caso:

>> test
time for (a''*(a*b)'')'': 1.717s
time for a*b*a: 1.7153s
error = 1.0914e-11

Me gustaría saber qué está causando este comportamiento y cómo realizar a*b*a o (a''*(a*b)'')'' o ambos en el menor tiempo (es decir, ~ 0.4s en lugar de ~ 1.7s) dentro de una función matlab en lugar de dentro de un script.

Answer 1

Me puse en contacto con el soporte técnico de Mathworks y Rylan finalmente arrojó algo de luz sobre este tema. (¡Gracias Rylan!) Su respuesta completa está abajo. El problema de la función frente a la secuencia de comandos parece estar relacionado con ciertas optimizaciones. Matlab se aplica automáticamente a las funciones (pero no a las secuencias de comandos) que no funcionan como se esperaba.

La respuesta de Rylan:

Gracias por su paciencia en este tema. He consultado con los desarrolladores de computación de la GPU de MATLAB para entender esto mejor.

Este problema se debe a optimizaciones internas realizadas por MATLAB cuando se encuentran algunas operaciones específicas como la multiplicación de matriz-matriz y la transposición. Algunas de estas optimizaciones pueden habilitarse específicamente cuando se ejecuta una función MATLAB (o una función anónima) en lugar de un script.

Cuando su código inicial se ejecutaba desde un script, no se realiza una optimización de transposición de matriz particular, lo que hace que la expresión ''res2'' sea más rápida que la expresión ''res1'':

  n = 2000;
  a=gpuArray(sprand(n,n,0.01)); 
  b=gpuArray(rand(n));

  tic;res1=a*b*a;wait(gpuDevice);toc                                         % Elapsed time is 0.884099 seconds.
  tic;res2=transpose(transpose(a)*transpose(a*b));wait(gpuDevice);toc        % Elapsed time is 0.068855 seconds.

Sin embargo, cuando el código anterior se coloca en un archivo de función MATLAB, se realiza una optimización adicional de los tiempos de transposición matricial que hace que la expresión ''res2'' pase por una ruta de código diferente (y una llamada a la función de biblioteca CUDA diferente) en comparación con la misma línea Llamado desde un script. Por lo tanto, esta optimización genera resultados más lentos para la línea ''res2'' cuando se llama desde un archivo de función.

Para evitar que ocurra este problema en un archivo de función, las operaciones de transposición y multiplicación deberían dividirse de una manera que impida que MATLAB aplique esta optimización. Separar cada cláusula dentro de la declaración ''res2'' parece ser suficiente para esto:

  tic;i1=transpose(a);i2=transpose(a*b);res3=transpose(i1*i2);wait(gpuDevice);toc      % Elapsed time is 0.066446 seconds.

En la línea anterior, ''res3'' se genera a partir de dos matrices intermedias: ''i1'' e ''i2''. El rendimiento (en mi sistema) parece estar a la par con el de la expresión ''res2'' cuando se ejecuta desde un script; Además, la expresión ''res3'' también muestra un rendimiento similar cuando se ejecuta desde un archivo de función MATLAB. Sin embargo, tenga en cuenta que se puede utilizar memoria adicional para almacenar la copia transpuesta de la matriz inicial. Déjeme saber si ve un comportamiento de rendimiento diferente en su sistema, y ​​puedo investigar esto más a fondo.

Además, la operación ''res3'' muestra un rendimiento más rápido cuando se mide con la función ''gputimeit'' también. Consulte el archivo ''testscript2.m'' adjunto para obtener más información sobre esto. También he adjuntado ''test_v2.m'', que es una modificación de la función ''test.m'' en tu publicación de desbordamiento de pila.

Gracias por informar de este problema a mí. Me gustaría disculparme por cualquier inconveniente causado por este problema. He creado un informe de error interno para notificar a los desarrolladores de MATLAB sobre este comportamiento. Pueden proporcionar una solución para esto en una futura versión de MATLAB.

Ya que tenía una pregunta adicional sobre la comparación del rendimiento del código de GPU usando ''gputimeit'' en lugar de ''tic'' y ''toc'', solo quería ofrecer una sugerencia que los desarrolladores de informática de MATLAB GPU habían mencionado anteriormente. En general, también es bueno llamar a ''esperar (gpuDevice)'' antes de las declaraciones ''tic'' para garantizar que las operaciones de GPU de las líneas anteriores no se superpongan en la medición de la siguiente línea. Por ejemplo, en las siguientes líneas:

  b=gpuArray(rand(n));
  tic; res1=a*b*a; wait(gpuDevice); toc  

Si no se llama el ''wait (gpuDevice)'' antes del ''tic'', parte del tiempo que se tarda en construir la matriz ''b'' desde la línea anterior puede superponerse y contarse en el tiempo necesario para ejecutar la expresión ''res1''. Esto sería preferible en su lugar:

  b=gpuArray(rand(n));
  wait(gpuDevice); tic; res1=a*b*a; wait(gpuDevice); toc  

Aparte de esto, no veo ningún problema específico en la forma en que está utilizando las funciones ''tic'' y ''toc''. Sin embargo, tenga en cuenta que el uso de ''gputimeit'' generalmente se recomienda sobre el uso de ''tic'' y ''toc'' directamente para el perfilado relacionado con la GPU.

Seguiré adelante y cerraré este caso por ahora, pero hágame saber si tiene más preguntas al respecto.

%testscript2.m
n = 2000;
a = gpuArray(sprand(n, n, 0.01)); 
b = gpuArray(rand(n)); 

gputimeit(@()transpose_mult_fun(a, b))
gputimeit(@()transpose_mult_fun_2(a, b))

function out = transpose_mult_fun(in1, in2)

i1 = transpose(in1);
i2 = transpose(in1*in2);

out = transpose(i1*i2);

end

function out = transpose_mult_fun_2(in1, in2)

out = transpose(transpose(in1)*transpose(in1*in2));

end

.

function test_v2

clear

%% transposed expression
n = 2000;
rng(''default'');rng(1);
a = sprand(n, n, 0.1);
b = rand(n, n);
a = gpuArray(a);
b = gpuArray(b);

tic;
c1 = gather(transpose( transpose(a) * transpose(a * b) ));

disp([''time for (a''''*(a*b)'''')'''': '' , num2str(toc),''s''])

clearvars -except c1

%% non-transposed expression
rng(''default'');
rng(1)
n = 2000;
a = sprand(n, n, 0.1);
b = rand(n, n);
a = gpuArray(a);
b = gpuArray(b);

tic;
c2 = gather(a * b * a);

disp([''time for a*b*a: '' , num2str(toc),''s''])
disp([''error = '',num2str(max(max(abs(c1-c2))))])

%% sliced equivalent
rng(''default'');
rng(1)
n = 2000;
a = sprand(n, n, 0.1);
b = rand(n, n);
a = gpuArray(a);
b = gpuArray(b);

tic;
intermediate1 = transpose(a);
intermediate2 = transpose(a * b);
c3 = gather(transpose( intermediate1 * intermediate2 ));

disp([''time for split equivalent: '' , num2str(toc),''s''])
disp([''error = '',num2str(max(max(abs(c1-c3))))])

end

Answer 2

Parece que hay un problema con la multiplicación de dos matrices dispersas en GPU. el tiempo para dispersar por matriz completa es más de 1000 veces más rápido que dispersado por dispersos. Un ejemplo simple:

str={''sparse*sparse'',''sparse*full''};
for ii=1:2
    rng(1);
    a=sprand(3000,3000,0.1);
    b=sprand(3000,3000,0.1);
    if ii==2
        b=full(b);
    end
    a=gpuArray(a);
    b=gpuArray(b);
    tic
    c=a*b;
    disp([''time for '',str{ii},'': '' , num2str(toc),''s''])
end

En tu contexto, es la última multiplicación la que lo hace. para demostrar que sustituyo a con un duplicado c , y multiplico por él dos veces, una vez como dispersa y una vez como matriz completa.

str={''a*b*a'',''a*b*full(a)''};
for ii=1:2
    %rng(''default'');
    rng(1)
    a=sprand(3000,3000,0.1);
    b=rand(3000,3000);
    rng(1)
    c=sprand(3000,3000,0.1);
    if ii==2
        c=full(c);
    end
    a=gpuArray(a);
    b=gpuArray(b);
    c=gpuArray(c);
    tic;
    c1{ii}=a*b*c;
    disp([''time for '',str{ii},'': '' , num2str(toc),''s''])
end
disp([''error = '',num2str(max(max(abs(c1{1}-c1{2}))))])

Puede que me equivoque, pero mi conclusión es que a * b * a implica la multiplicación de dos matrices dispersas (aya de nuevo) y no se trata bien, mientras que el enfoque de transposición () divide el proceso en dos etapas, en ninguna de que hay dos matrices dispersas.

Answer 3

EDIT 2 Podría haber tenido razón, vea esta otra respuesta

EDITAR : Utilizan MAGMA , que es la columna mayor. Mi respuesta no es válida, sin embargo, la dejaré aquí por un tiempo en caso de que pueda ayudar a resolver este extraño comportamiento.

La siguiente respuesta es incorrecta

Esta es mi suposición, no puedo decirle al 100% sin saber el código debajo de la capucha de MATLAB.

Hipótesis: el código de computación paralela de MATLAB utiliza las bibliotecas CUDA, no las propias.

Información importante

• MATLAB es columna principal y CUDA es fila mayor.
• No hay cosas como matrices 2D, solo matrices 1D con 2 índices

¿Por qué importa esto? Bien, porque CUDA es un código altamente optimizado que usa la estructura de memoria para maximizar los aciertos de caché por kernel (la operación más lenta en las GPU es la memoria de lectura). Esto significa que un código de multiplicación de matriz CUDA estándar explotará el orden de las lecturas de la memoria para asegurarse de que sean adyacentes. Sin embargo, lo que es la memoria adyacente en la fila principal no está en la columna principal.

Por lo tanto, hay 2 soluciones a esto como alguien que escribe software

1. Escribe tus propias bibliotecas de álgebra de columna principal en CUDA
2. Tome cada entrada / salida de MATLAB y transpórtela (es decir, convierta de columna principal a fila mayor)

Han hecho el punto 2, y suponiendo que hay un compilador JIT inteligente para MATLAB, caja de herramientas de procesamiento paralelo (suposición razonable), para el segundo caso, toma a y b , las transpone, realiza los cálculos y transpone la salida cuando se gather .

Sin embargo, en el primer caso ya no es necesario transponer la salida, ya que ya está transpuesta internamente y el JIT lo capta, por lo que, en lugar de llamar a gather(transpose( XX )) , simplemente omite la transposición de salida. Lo mismo con la transpose(a*b) . Tenga en cuenta que transpose(a*b)=transpose(b)*transpose(a) , por lo que de repente no se necesitan transposiciones (todas se omiten internamente). Una transposición es una operación costosa.

De hecho, hay algo extraño aquí: convertir el código en una función de repente lo hace lento. Mi mejor conjetura es que debido a que el JIT se comporta de manera diferente en diferentes situaciones, no atrapa todas estas cosas de transpose interior y simplemente realiza todas las operaciones, perdiendo la velocidad.

Observación interesante: se necesita el mismo tiempo en la CPU que en la GPU para hacer a*b*a en mi PC.