O arrayfun pode ser significativamente mais lento que um loop explícito no matlab. Por quê?

Considere o seguinte teste de velocidade simples 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 

Na minha máquina (Matlab 2011b no Linux Mint 12), a saída deste teste é:

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

O que?!? arrayfun , embora seja reconhecidamente uma solução mais limpa, é uma ordem de grandeza mais lenta. O que está acontecendo aqui?

Além disso, fiz um estilo semelhante de teste para o cellfun e descobri que ele é cerca de 3 vezes mais lento que um loop explícito. Mais uma vez, esse resultado é o oposto do que eu esperava.

Minha pergunta é: Por que o arrayfun e o cellfun são muito mais lentos? E, dado isso, há alguma boa razão para usá-los (além de fazer o código parecer bom)?

Nota: Eu estou falando sobre a versão padrão do arrayfun aqui, NÃO a versão GPU da checkbox de ferramentas de parallel processing.

EDIT: Só para ficar claro, estou ciente de que Func1 acima pode ser vetorizado como apontado por Oli. Eu escolhi apenas porque ele produz um teste de velocidade simples para os propósitos da questão real.

EDIT: Após a sugestão de grungetta, eu refiz o teste com feature accel off . Os resultados são:

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

Em outras palavras, parece que uma grande parte da diferença é que o acelerador JIT faz um trabalho muito melhor de acelerar o loop for explícito do que o arrayfun . Isso parece estranho para mim, pois o arrayfun realmente fornece mais informações, ou seja, seu uso revela que a ordem das chamadas para Func1 não importa. Além disso, notei que se o acelerador JIT é ligado ou desligado, meu sistema usa apenas uma CPU …

Você pode obter a ideia executando outras versões do seu código. Considere explicitamente escrever os cálculos, em vez de usar uma function em seu loop

 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 

Hora de calcular no meu computador:

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

Agora, enquanto a solução totalmente ‘vetorizada’ é claramente a mais rápida, você pode ver que definir uma function para ser chamada para cada input x é uma sobrecarga enorme . Apenas explicitamente escrever o cálculo nos deu aceleração de fator 5. Eu acho que isso mostra que o compilador JIT do MATLAB não suporta funções embutidas . De acordo com a resposta do gnovice, é melhor escrever uma function normal do que uma anônima. Tente.

Próximo passo – remova (vetorize) o laço 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. 

Outro fator 5 aceleração: há algo nessas declarações dizendo que você deve evitar loops no MATLAB … Ou existe mesmo? Dê uma olhada nisso então

 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. 

Muito mais perto da versão ‘totalmente’ vectorizada. O Matlab armazena matrizes em coluna. Você deve sempre (quando possível) estruturar seus cálculos para serem vetorizados ’em coluna’.

Nós podemos voltar para o Soln3 agora. A ordem de loop é ‘row-wise’. Vamos mudar isso

 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. 

Melhor, mas ainda é muito ruim. Único loop – bom. Laço duplo – ruim. Eu acho que o MATLAB fez algum trabalho decente em melhorar o desempenho dos loops, mas ainda assim a sobrecarga do loop está lá. Se você tivesse algum trabalho pesado dentro, você não notaria. Mas como esse cálculo é limitado pela largura de banda da memory, você vê a sobrecarga do loop. E você vai ver ainda mais claramente a sobrecarga de chamar Func1 lá.

Então, o que há com o arrayfun? Não há nenhuma function interna ali, então muita sobrecarga. Mas por que muito pior do que um loop duplo nested? Na verdade, o tópico de usar o cellfun / arrayfun foi amplamente discutido várias vezes (por exemplo, aqui , aqui , aqui e aqui ). Essas funções são simplesmente lentas, você não pode usá-las para cálculos finos. Você pode usá-los para brevidade de código e conversões extravagantes entre células e matrizes. Mas a function precisa ser mais pesada do que você escreveu:

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

Note que Soln7 é uma célula agora .. às vezes isso é útil. O desempenho do código é muito bom agora, e se você precisar de célula como saída, não precisará converter sua matriz depois de ter usado a solução totalmente vetorizada.

Então, por que é arrayfun mais lento que uma estrutura de loop simples? Infelizmente, é impossível dizermos com certeza, já que não há código-fonte disponível. Você só pode imaginar que, como arrayfun é uma function de propósito geral, que lida com todos os tipos de estruturas de dados e argumentos diferentes, não é necessariamente muito rápido em casos simples, que você pode expressar diretamente como ninhos de loop. De onde vem a sobrecarga não podemos saber. A sobrecarga poderia ser evitada por uma implementação melhor? Talvez não. Mas, infelizmente, a única coisa que podemos fazer é estudar o desempenho para identificar os casos, nos quais ele funciona bem e aqueles em que não funciona.

Atualização Como o tempo de execução deste teste é curto, para obter resultados confiáveis, adicionei agora um loop em torno dos testes:

 for i=1:1000 % compute end 

Algumas vezes abaixo:

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

Você vê que o arrayfun ainda é ruim, mas pelo menos três ordens de magnitude pior que a solução vetorizada. Por outro lado, um único loop com cálculos de coluna é tão rápido quanto a versão totalmente vetorizada … Isso foi feito em uma única CPU. Resultados para Soln5 e Soln7 não mudam se eu mudar para 2 núcleos – No Soln5 eu teria que usar um parfor para obtê-lo em paralelo. Esqueça o speedup … O Soln7 não roda em paralelo porque o arrayfun não roda em paralelo. A versão vetorizada de Olis, por outro lado:

 Oli 5.508085 seconds. 

Isto porque!!!!

 x = randn(T, N); 

não é tipo gpuarray ;

Tudo que você precisa fazer é

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