arrayfun puede ser significativamente más lento que un bucle 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 , si bien es cierto que es una solución de aspecto más limpio, es un orden de magnitud más lento. ¿Que esta pasando aqui?

Además, hice un estilo de prueba similar para cellfun y encontré que es aproximadamente 3 veces más lento que un loop explícito. De nuevo, este resultado es el opuesto de lo que esperaba.

Mi pregunta es: ¿por qué arrayfun y cellfun son mucho más lentos? Y dado esto, ¿hay alguna buena razón para usarlos (aparte de hacer que el código se vea bien)?

Nota: Estoy hablando de la versión estándar de arrayfun aquí, NO de la versión de GPU de la caja de herramientas de parallel processing.

EDITAR: Para ser claros, soy consciente de que Func1 arriba puede ser vectorizado como lo señala Oli. Solo lo elegí porque arroja 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 la feature accel off . Los resultados son:

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

En otras palabras, parece que una gran parte de la diferencia es que el acelerador de JIT hace un trabajo mucho mejor al acelerar el ciclo for explícito for lo que hace arrayfun . Esto me parece extraño, ya que arrayfun proporciona más información, es decir, su uso revela que el orden de las llamadas a Func1 no importa. Además, noté que si el acelerador JIT está encendido o apagado, mi sistema solo usa una CPU …

Puede obtener la 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 

Tiempo para 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 la definición de una función que debe invocarse para cada entrada x es una gran sobrecarga. El hecho de escribir explícitamente el cálculo nos permitió acelerar el factor 5. Supongo que esto muestra que el comstackdor 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 anónima. Intentalo.

Siguiente paso – eliminar (vectorizar) el ciclo 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 acelerador de factor 5: hay algo en esas afirmaciones que dice que debes evitar los bucles en MATLAB … ¿O realmente hay? Eche 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 ‘totalmente’ vectorizada. Matlab almacena las matrices en columnas. Siempre debe (cuando sea posible) estructurar sus cálculos para que sean vectorizados ‘en columnas’.

Podemos volver a Soln3 ahora. El orden de bucle es ‘en hilera’. 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 aún muy malo. Bucle único – bien. Doble bucle – malo. Supongo que MATLAB hizo un trabajo decente para mejorar el rendimiento de los bucles, pero aún así el ciclo de sobrecarga está ahí. Si tuviera un trabajo más pesado adentro, no lo notaría. Pero dado que este cálculo está limitado al ancho de banda de la memoria, sí se ve el ciclo por encima. Y verá aún más claramente la sobrecarga de llamar a Func1 allí.

Entonces, ¿qué pasa con arrayfun? Ninguna función entra allí, así que mucha sobrecarga. Pero, ¿por qué tanto peor que un bucle nested doble? En realidad, el tema del uso de cellfun / arrayfun ha sido extensamente 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 necesita ser más pesada que lo que escribiste:

 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 Soln7 es una celda ahora … a veces eso es útil. El rendimiento del código es bastante bueno ahora, y si necesita una celda como salida, no necesita convertir su matriz después de haber utilizado la solución completamente vectorizada.

Entonces, ¿por qué arrayfun es más lento que una estructura simple de bucle? Desafortunadamente, es imposible para nosotros decirlo con seguridad, ya que no hay un código fuente disponible. Solo se 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ápido en casos simples, que puede express directamente como bucle anida. ¿De dónde viene la sobrecarga que no podemos saber? ¿Se podría evitar la sobrecarga mediante 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 aquellos, donde 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 dado a continuación:

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

Usted ve 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 totalmente vectorizada … Todo 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 hacerlo paralelizar. Olvídese de la aceleración … Soln7 no se ejecuta en paralelo porque arrayfun no se ejecuta en paralelo. La versión vectorizada de Olis por otro lado:

 Oli 5.508085 seconds. 

Eso porque !!!!

 x = randn(T, N); 

no es tipo gpuarray ;

Todo lo que necesitas hacer es

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