arrayfun pode ser significativamente mais lento do 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), o resultado deste teste é:

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

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

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

Minha pergunta é: Por que são arrayfune cellfuntão mais lentos? E, considerando isso, há algum bom motivo para usá-los (além de fazer o código parecer bom)?

Observação: estou falando sobre a versão padrão arrayfunaqui, NÃO a versão da GPU da caixa de ferramentas de processamento paralelo.

EDIT: Só para ficar claro, estou ciente de que Func1acima pode ser vetorizado conforme apontado por Oli. Eu só o escolhi porque ele produz um teste de velocidade simples para os propósitos da pergunta real.

EDIT: Seguindo a sugestão de grungetta, 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, pareceria que uma grande parte da diferença é que o acelerador JIT faz um trabalho muito melhor em acelerar o forloop explícito do que faz arrayfun. Isso me parece estranho, pois arrayfunna verdade fornece mais informações, ou seja, seu uso revela que a ordem das chamadas para Func1não importa. Além disso, observei que, quer o acelerador JIT esteja ativado ou desativado, meu sistema usa apenas uma CPU ...

Você pode ter uma ideia executando outras versões do seu código. Considere escrever explicitamente os cálculos, em vez de usar uma função 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, embora a solução totalmente 'vetorizada' seja claramente a mais rápida, você pode ver que definir uma função a ser chamada para cada entrada x é uma grande sobrecarga. O simples fato de escrever explicitamente o cálculo nos fez aumentar a velocidade do fator 5. Eu acho que isso mostra que o compilador MATLABs JIT não oferece suporte a funções embutidas . De acordo com a resposta de gnovice lá, é realmente melhor escrever uma função normal do que uma anônima. Tente.

Próxima etapa - remover (vetorizar) o loop 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 de aceleração: há algo nessas declarações dizendo que você deve evitar loops no MATLAB ... Ou será que há 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 próximo da versão 'totalmente' vetorizada. Matlab armazena matrizes em colunas. Você deve sempre (quando possível) estruturar seus cálculos para serem vetorizados 'em colunas'.

Podemos voltar para Soln3 agora. A ordem do loop é 'por linha'. 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. Loop único - bom. Loop duplo - ruim. Eu acho que o MATLAB fez um trabalho decente para melhorar o desempenho dos loops, mas ainda assim a sobrecarga do loop está lá. Se você tivesse algum trabalho mais pesado por dentro, não notaria. Mas como esse cálculo é limitado pela largura de banda da memória, você vê a sobrecarga do loop. E você verá ainda mais claramente a sobrecarga de chamar Func1 lá.

Então, o que há com arrayfun? Nenhuma função inlinig lá também, então muita sobrecarga. Mas por que é tão pior do que um loop aninhado duplo? Na verdade, o tema do uso cellfun / arrayfun foi extensamente discutido muitas vezes (por exemplo, aqui , aqui , aqui e aqui ). Essas funções são simplesmente lentas, você não pode usá-las para cálculos de granulação fina. Você pode usá-los para abreviar o código e conversões sofisticadas entre células e matrizes. Mas a função precisa ser mais pesada do que o que você escreveu:

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

Soln7  0.016786 seconds.

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

Então, por que arrayfun é mais lento do que uma estrutura de loop simples? Infelizmente, é impossível para nós ter certeza, uma vez que não há código-fonte disponível. Você só pode supor que, como arrayfun é uma função de propósito geral, que lida com todos os tipos de estruturas de dados e argumentos diferentes, ela não é necessariamente muito rápida 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 em que 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 não três ordens de magnitude pior do que a solução vetorial. Por outro lado, um único loop com cálculos em colunas é tão rápido quanto a versão totalmente vetorizada ... Tudo isso foi feito em uma única CPU. Os resultados para Soln5 e Soln7 não mudam se eu mudar para 2 núcleos - Em Soln5, eu teria que usar um parfor para colocá-lo em paralelo. Esqueça o speedup ... Soln7 não roda em paralelo porque arrayfun não roda em paralelo. Versão vetorizada Olis por outro lado:

Oli  5.508085 seconds.

Isto porque!!!!

x = randn(T, N); 

não é gpuarraytipo;

Tudo que você precisa fazer é

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