arrayfun может быть значительно медленнее, чем явный цикл в Matlab. Зачем?

Рассмотрим следующий простой тест скорости 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

На моей машине (Matlab 2011b на Linux Mint 12) результат этого теста:

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

Что за?!? arrayfunхотя и выглядит более чистым, но на порядок медленнее. Что здесь происходит?

Кроме того, я проделал похожий стиль тестирования cellfunи обнаружил, что он примерно в 3 раза медленнее, чем явный цикл. Опять же, этот результат противоположен тому, что я ожидал.

Мой вопрос: Почему arrayfunи cellfunтак много медленнее? И, учитывая это, есть ли какие-нибудь веские причины использовать их (кроме как для того, чтобы код выглядел хорошо)?

Примечание: здесь я говорю о стандартной версии arrayfun, а НЕ о версии для графического процессора из набора инструментов параллельной обработки.

РЕДАКТИРОВАТЬ: Чтобы быть ясным, я знаю, что Func1выше можно векторизовать, как указано Оли. Я выбрал его только потому, что он дает простой тест скорости для целей фактического вопроса.

РЕДАКТИРОВАТЬ: Следуя предложению grungetta, я повторно провел тест с feature accel off. Результат:

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

Другими словами, может показаться, что большая часть разницы заключается в том, что ускоритель JIT выполняет гораздо лучшую работу по ускорению явного forцикла, чем он arrayfun. Мне это кажется странным, поскольку на arrayfunсамом деле предоставляет больше информации, т. Е. Его использование показывает, что порядок вызовов Func1не имеет значения. Кроме того, я заметил, что независимо от того, включен ли JIT-ускоритель или выключен, моя система всегда использует только один процессор ...

Вы можете получить представление, запустив другие версии вашего кода. Рассмотрите возможность явной записи вычислений вместо использования функции в вашем цикле

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

Время вычислить на моем компьютере:

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

Теперь, хотя полностью «векторизованное» решение явно является самым быстрым, вы можете видеть, что определение функции, которая будет вызываться для каждой записи x, - это огромные накладные расходы. Просто явная запись вычислений дала нам ускорение в 5 раз. Я думаю, это показывает, что JIT-компилятор MATLABs не поддерживает встроенные функции . Согласно ответу gnovice там на самом деле лучше написать нормальную функцию, чем анонимную. Попытайся.

Следующий шаг - удалить (векторизовать) внутренний цикл:

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.

Еще один фактор 5 ускорения: в этих утверждениях есть что-то, что вам следует избегать циклов в MATLAB ... Или это действительно так? Посмотри на это тогда

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.

Намного ближе к «полностью» векторизованной версии. Matlab хранит матрицы по столбцам. Вы всегда должны (по возможности) структурировать свои вычисления так, чтобы они были векторизованы «по столбцам».

Теперь мы можем вернуться к Soln3. Порядок циклов там «построчный». Давай изменим это

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.

Лучше, но все равно очень плохо. Одиночная петля - хорошо. Двойная петля - плохо. Я предполагаю, что MATLAB проделал приличную работу по повышению производительности циклов, но накладные расходы на цикл все же присутствуют. Если бы у вас внутри была более тяжелая работа, вы бы этого не заметили. Но поскольку это вычисление ограничено пропускной способностью памяти, вы видите накладные расходы цикла. И вы будете еще более ясно увидеть накладные расходы вызова FUNC1 там.

Так что случилось с arrayfun? Там тоже нет функции inlinig, поэтому много накладных расходов. Но почему это намного хуже, чем двойной вложенный цикл? На самом деле, тема использования cellfun / arrayfun широко обсуждалась много раз (например, здесь , здесь , здесь и здесь ). Эти функции просто медленные, вы не можете использовать их для таких мелкозернистых вычислений. Вы можете использовать их для краткости кода и необычных преобразований между ячейками и массивами. Но функция должна быть тяжелее того, что вы написали:

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

Soln7  0.016786 seconds.

Обратите внимание, что Soln7 теперь является ячейкой ... иногда это полезно. Производительность кода сейчас неплохая, и если вам нужна ячейка в качестве вывода, вам не нужно преобразовывать матрицу после использования полностью векторизованного решения.

Так почему же arrayfun медленнее, чем простая циклическая структура? К сожалению, мы не можем сказать наверняка, поскольку исходный код недоступен. Вы можете только догадываться, что, поскольку arrayfun - это функция общего назначения, которая обрабатывает все виды различных структур данных и аргументов, она не обязательно очень быстрая в простых случаях, которые вы можете напрямую выразить как вложения циклов. Откуда берутся накладные расходы, мы не можем знать. Можно ли избежать накладных расходов за счет лучшей реализации? Может быть нет. Но, к сожалению, единственное, что мы можем сделать, - это изучить производительность, чтобы определить, в каких случаях она работает хорошо, а в каких - нет.

Обновление Поскольку время выполнения этого теста невелико, для получения надежных результатов я добавил цикл вокруг тестов:

for i=1:1000
   % compute
end

Несколько раз указано ниже:

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

Вы видите, что arrayfun по-прежнему плохой, но как минимум не на три порядка хуже, чем векторизованное решение. С другой стороны, один цикл с вычислениями по столбцам выполняется так же быстро, как и полностью векторизованная версия ... Все это было сделано на одном процессоре. Результаты для Soln5 и Soln7 не меняются, если я переключаюсь на 2 ядра - в Soln5 мне пришлось бы использовать parfor для его распараллеливания. Забудьте об ускорении ... Soln7 не работает параллельно, потому что arrayfun не работает параллельно. С другой стороны, векторизованная версия Olis:

Oli  5.508085 seconds.

Это потому что !!!!

x = randn(T, N); 

не gpuarrayтип;

Все, что вам нужно сделать, это

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