arrayfun puede ser significativamente más lento que un ciclo explícito en matlab. ¿Por qué?

Considere la siguiente prueba de velocidad simple para arrayfun:

T = 4000;
N = 500;
x = randn(T, N);
Func1 = @(a) (3*a^2 + 2*a - 1);

tic
Soln1 = ones(T, N);
for t = 1:T
    for n = 1:N
        Soln1(t, n) = Func1(x(t, n));
    end
end
toc

tic
Soln2 = arrayfun(Func1, x);
toc

En mi máquina (Matlab 2011b en Linux Mint 12), el resultado de esta prueba es:

Elapsed time is 1.020689 seconds.
Elapsed time is 9.248388 seconds.

¿¡¿Que?!? arrayfun, aunque es una solución de apariencia más limpia, es un orden de magnitud más lenta. ¿Que esta pasando aqui?

Además, hice un estilo similar de prueba cellfuny descubrí que era aproximadamente 3 veces más lento que un bucle explícito. Nuevamente, este resultado es el opuesto al que esperaba.

Mi pregunta es: ¿Por qué son arrayfuny cellfunmucho más lentos? Y dado esto, ¿hay alguna buena razón para usarlos (además de hacer que el código se vea bien)?

Nota: estoy hablando de la versión estándar de arrayfunaquí, NO de la versión de GPU de la caja de herramientas de procesamiento paralelo.

EDITAR: Para que quede claro, soy consciente de que lo Func1anterior se puede vectorizar como lo señaló Oli. Solo lo elegí porque produce una prueba de velocidad simple para los propósitos de la pregunta real.

EDITAR: Siguiendo la sugerencia de grungetta, volví a hacer la prueba con feature accel off. Los resultados son:

Elapsed time is 28.183422 seconds.
Elapsed time is 23.525251 seconds.

En otras palabras, parecería que una gran parte de la diferencia es que el acelerador JIT acelera mucho mejor el forciclo explícito que lo hace arrayfun. Esto me parece extraño, ya que en arrayfunrealidad brinda más información, es decir, su uso revela que el orden de las llamadas Func1no importa. Además, noté que si el acelerador JIT está encendido o apagado, mi sistema solo usa una CPU ...

Puede hacerse una idea ejecutando otras versiones de su código. Considere escribir explícitamente los cálculos, en lugar de usar una función en su ciclo

tic
Soln3 = ones(T, N);
for t = 1:T
    for n = 1:N
        Soln3(t, n) = 3*x(t, n)^2 + 2*x(t, n) - 1;
    end
end
toc

Es hora de calcular en mi computadora:

Soln1  1.158446 seconds.
Soln2  10.392475 seconds.
Soln3  0.239023 seconds.
Oli    0.010672 seconds.

Ahora, mientras que la solución completamente 'vectorizada' es claramente la más rápida, puede ver que definir una función que se llamará para cada entrada x es una sobrecarga enorme . Simplemente escribir explícitamente el cálculo nos dio una aceleración del factor 5. Supongo que esto muestra que el compilador JIT de MATLAB no admite funciones en línea . Según la respuesta de gnovice allí, en realidad es mejor escribir una función normal en lugar de una anónima. Intentalo.

Siguiente paso: eliminar (vectorizar) el bucle interno:

tic
Soln4 = ones(T, N);
for t = 1:T
    Soln4(t, :) = 3*x(t, :).^2 + 2*x(t, :) - 1;
end
toc

Soln4  0.053926 seconds.

Otro factor 5 de aceleración: hay algo en esas declaraciones que dice que debe evitar los bucles en MATLAB ... ¿O realmente? Echa un vistazo a esto entonces

tic
Soln5 = ones(T, N);
for n = 1:N
    Soln5(:, n) = 3*x(:, n).^2 + 2*x(:, n) - 1;
end
toc

Soln5   0.013875 seconds.

Mucho más cerca de la versión "completamente" vectorizada. Matlab almacena matrices por columnas. Siempre debe (cuando sea posible) estructurar sus cálculos para que sean vectorizados 'en columnas'.

Podemos volver a Soln3 ahora. El orden de los bucles es "por filas". Vamos a cambiarlo

tic
Soln6 = ones(T, N);
for n = 1:N
    for t = 1:T
        Soln6(t, n) = 3*x(t, n)^2 + 2*x(t, n) - 1;
    end
end
toc

Soln6  0.201661 seconds.

Mejor, pero muy malo. Bucle único - bueno. Doble bucle - malo. Supongo que MATLAB hizo un trabajo decente para mejorar el rendimiento de los bucles, pero aún así la sobrecarga del bucle está ahí. Si tuvieras un trabajo más pesado adentro, no lo notarías. Pero dado que este cálculo está limitado al ancho de banda de la memoria, se ve la sobrecarga del bucle. Y verá aún más claramente la sobrecarga de llamar a Func1 allí.

Entonces, ¿qué pasa con arrayfun? Allí tampoco hay ninguna función, por lo que hay muchos gastos generales. Pero, ¿por qué es mucho peor que un bucle anidado doble? En realidad, el tema de la utilización de cellfun / arrayfun ha sido ampliamente discutido muchas veces (por ejemplo, aquí , aquí , aquí y aquí ). Estas funciones son simplemente lentas, no puede usarlas para cálculos tan finos. Puede usarlos para la brevedad del código y las conversiones sofisticadas entre celdas y matrices. Pero la función debe ser más pesada de lo que escribió:

tic
Soln7 = arrayfun(@(a)(3*x(:,a).^2 + 2*x(:,a) - 1), 1:N, 'UniformOutput', false);
toc

Soln7  0.016786 seconds.

Tenga en cuenta que ahora Soln7 es una celda ... a veces eso es útil. El rendimiento del código es bastante bueno ahora, y si necesita celdas como salida, no necesita convertir su matriz después de haber usado la solución completamente vectorizada.

Entonces, ¿por qué arrayfun es más lento que una estructura de bucle simple? Desafortunadamente, es imposible para nosotros decirlo con certeza, ya que no hay un código fuente disponible. Solo puede adivinar que, dado que arrayfun es una función de propósito general, que maneja todo tipo de estructuras de datos y argumentos diferentes, no es necesariamente muy rápida en casos simples, que puede expresar directamente como nidos de bucle. No podemos saber de dónde viene la sobrecarga. ¿Se podrían evitar los gastos generales con una mejor implementación? Tal vez no. Pero, lamentablemente, lo único que podemos hacer es estudiar el rendimiento para identificar los casos en los que funciona bien y en los que no.

Actualización Dado que el tiempo de ejecución de esta prueba es corto, para obtener resultados confiables agregué ahora un ciclo alrededor de las pruebas:

for i=1:1000
   % compute
end

Algunas veces se dan a continuación:

Soln5   8.192912 seconds.
Soln7  13.419675 seconds.
Oli     8.089113 seconds.

Verá que el arrayfun sigue siendo malo, pero al menos no tres órdenes de magnitud peor que la solución vectorizada. Por otro lado, un solo bucle con cálculos en columnas es tan rápido como la versión completamente vectorizada ... Todo eso se hizo en una sola CPU. Los resultados para Soln5 y Soln7 no cambian si cambio a 2 núcleos - En Soln5 tendría que usar un parfor para paralelizarlo. Olvídese de la aceleración ... Soln7 no se ejecuta en paralelo porque arrayfun no se ejecuta en paralelo. Versión vectorizada de Olis por otro lado:

Oli  5.508085 seconds.

Eso es porque !!!!

x = randn(T, N); 

no es gpuarraytipo;

Todo lo que necesitas hacer es

x = randn(T, N,'gpuArray');