Ответ 1
Это нормальное и ожидаемое поведение. Это слишком упрощено для применения инструкции "избегать для циклов с numpy" . Если вы имеете дело с внутренними петлями, то это почти всегда верно. Но в случае внешних петель (например, в вашем случае) есть гораздо больше исключений. Особенно, если альтернативой является использование широковещательной передачи, поскольку это ускоряет работу, используя намного больше.
Просто добавьте немного фона для этого "избегать циклов с numpy" :
Массивы NumPy сохраняются в виде непрерывных массивов с c. Python int не совпадает с C int! Поэтому всякий раз, когда вы перебираете каждый элемент в массиве, вам нужно подключить элемент из массива, преобразовать его в Python int, а затем делать все, что вы хотите с ним сделать, и, наконец, вам может понадобиться снова преобразовать его в переменную integer (называемый бокс и распаковка значения). Например, вы хотите sum элементы в массиве с помощью Python:
import numpy as np
arr = np.arange(1000)
%%timeit
acc = 0
for item in arr:
acc += item
# 1000 loops, best of 3: 478 µs per loop
Лучше использовать numpy:
%timeit np.sum(arr)
# 10000 loops, best of 3: 24.2 µs per loop
Даже если вы нажмете цикл в код Python C, вы далеко от производительности numpy:
%timeit sum(arr)
# 1000 loops, best of 3: 387 µs per loop
Из этого правила могут быть исключения, но они будут очень скудными. По крайней мере, пока существует некоторая эквивалентная функциональность numpy. Поэтому, если вы будете перебирать отдельные элементы, вы должны использовать numpy.
Иногда достаточно простого цикла python. Это не слишком рекламируемые, но функции numpy имеют огромные накладные расходы по сравнению с функциями Python. Например, рассмотрим 3-элементный массив:
arr = np.arange(3)
%timeit np.sum(arr)
%timeit sum(arr)
Какой из них будет быстрее?
Решение: функция Python работает лучше, чем решение numpy:
# 10000 loops, best of 3: 21.9 µs per loop <- numpy
# 100000 loops, best of 3: 6.27 µs per loop <- python
Но что это имеет отношение к вашему примеру? На самом деле это не так, потому что вы всегда используете numpy-функции для массивов (не отдельные элементы и даже не несколько элементов), поэтому ваш внутренний цикл уже использует оптимизированные функции. Вот почему обе выполняют примерно одинаковые (+/- коэффициент 10 с очень небольшим количеством элементов до коэффициента 2 примерно с 500 элементами). Но на самом деле это не слишком сложный цикл, это вызов функции вызова!
Решение вашей петли
Используя line-profiler и resolution = 100:
def fun_func(tim, prec, values):
for i, ti in enumerate(tim):
values[i] = np.sum(np.sin(prec * ti))
%lprun -f fun_func fun_func(tim, prec, values)
Line # Hits Time Per Hit % Time Line Contents
==============================================================
1 def fun_func(tim, prec, values):
2 101 752 7.4 5.7 for i, ti in enumerate(tim):
3 100 12449 124.5 94.3 values[i] = np.sum(np.sin(prec * ti))
95% проводится внутри цикла, я даже разбил тело цикла на несколько частей, чтобы проверить это:
def fun_func(tim, prec, values):
for i, ti in enumerate(tim):
x = prec * ti
x = np.sin(x)
x = np.sum(x)
values[i] = x
%lprun -f fun_func fun_func(tim, prec, values)
Line # Hits Time Per Hit % Time Line Contents
==============================================================
1 def fun_func(tim, prec, values):
2 101 609 6.0 3.5 for i, ti in enumerate(tim):
3 100 4521 45.2 26.3 x = prec * ti
4 100 4646 46.5 27.0 x = np.sin(x)
5 100 6731 67.3 39.1 x = np.sum(x)
6 100 714 7.1 4.1 values[i] = x
Потребители времени np.multiply, np.sin, np.sum здесь, как вы можете легко проверить, сравнивая их время на вызов с их накладными расходами:
arr = np.ones(1, float)
%timeit np.sum(arr)
# 10000 loops, best of 3: 22.6 µs per loop
Итак, как только накладные расходы на коллатульную функцию малы по сравнению с временем выполнения расчета, вы будете иметь схожие режимы работы. Даже со 100 пунктами вы достаточно близки к накладным расходам. Трюк заключается в том, чтобы знать, в какой момент они безупречны. При использовании 1000 элементов накладные расходы по-прежнему значительны:
%lprun -f fun_func fun_func(tim, prec, values)
Line # Hits Time Per Hit % Time Line Contents
==============================================================
1 def fun_func(tim, prec, values):
2 1001 5864 5.9 2.4 for i, ti in enumerate(tim):
3 1000 42817 42.8 17.2 x = prec * ti
4 1000 119327 119.3 48.0 x = np.sin(x)
5 1000 73313 73.3 29.5 x = np.sum(x)
6 1000 7287 7.3 2.9 values[i] = x
Но с resolution = 5000 накладные расходы довольно низки по сравнению с временем выполнения:
Line # Hits Time Per Hit % Time Line Contents
==============================================================
1 def fun_func(tim, prec, values):
2 5001 29412 5.9 0.9 for i, ti in enumerate(tim):
3 5000 388827 77.8 11.6 x = prec * ti
4 5000 2442460 488.5 73.2 x = np.sin(x)
5 5000 441337 88.3 13.2 x = np.sum(x)
6 5000 36187 7.2 1.1 values[i] = x
Когда вы потратили 500 долларов в каждом звонке np.sin, вам больше не нужно беспокоиться о 20us накладных расходах.
Вероятнее всего, должно быть слово предостережения: line_profiler, вероятно, включает в себя некоторые дополнительные накладные расходы на строку, возможно, и на вызов функции, поэтому точка, в которой служебные вызовы функции становятся небрежными, может быть ниже.
Ваше широковещательное решение
Я начал с профилирования первого решения, сделайте то же самое со вторым решением:
def fun_func(tim, prec, values):
x = tim[:, np.newaxis]
x = x * prec
x = np.sin(x)
x = np.sum(x, axis=1)
return x
Снова с помощью line_profiler с resolution=100:
%lprun -f fun_func fun_func(tim, prec, values)
Line # Hits Time Per Hit % Time Line Contents
==============================================================
1 def fun_func(tim, prec, values):
2 1 27 27.0 0.5 x = tim[:, np.newaxis]
3 1 638 638.0 12.9 x = x * prec
4 1 3963 3963.0 79.9 x = np.sin(x)
5 1 326 326.0 6.6 x = np.sum(x, axis=1)
6 1 4 4.0 0.1 return x
Это уже значительно превышает накладное время, и, таким образом, мы оказываемся в 10 раз быстрее по сравнению с циклом.
Я также сделал профилирование для resolution=1000:
Line # Hits Time Per Hit % Time Line Contents
==============================================================
1 def fun_func(tim, prec, values):
2 1 28 28.0 0.0 x = tim[:, np.newaxis]
3 1 17716 17716.0 14.6 x = x * prec
4 1 91174 91174.0 75.3 x = np.sin(x)
5 1 12140 12140.0 10.0 x = np.sum(x, axis=1)
6 1 10 10.0 0.0 return x
и precision=5000:
Line # Hits Time Per Hit % Time Line Contents
==============================================================
1 def fun_func(tim, prec, values):
2 1 34 34.0 0.0 x = tim[:, np.newaxis]
3 1 333685 333685.0 11.1 x = x * prec
4 1 2391812 2391812.0 79.6 x = np.sin(x)
5 1 280832 280832.0 9.3 x = np.sum(x, axis=1)
6 1 14 14.0 0.0 return x
Размер 1000 по-прежнему выполняется быстрее, но, как мы видели, в решении loop все еще не было пренебрежимо неверно. Но для resolution = 5000 время, затрачиваемое на каждом шаге, почти идентично (некоторые из них немного медленнее, другие быстрее, но в целом довольно похожи)
Другой эффект заключается в том, что фактическое вещание, когда вы выполняете умножение, становится значительным. Даже при использовании очень умных решений numpy это все еще включает в себя некоторые дополнительные вычисления. Для resolution=10000 вы видите, что трансляция вещания начинает занимать больше "% времени" по отношению к решению цикла:
Line # Hits Time Per Hit % Time Line Contents
==============================================================
1 def broadcast_solution(tim, prec, values):
2 1 37 37.0 0.0 x = tim[:, np.newaxis]
3 1 1783345 1783345.0 13.9 x = x * prec
4 1 9879333 9879333.0 77.1 x = np.sin(x)
5 1 1153789 1153789.0 9.0 x = np.sum(x, axis=1)
6 1 11 11.0 0.0 return x
Line # Hits Time Per Hit % Time Line Contents
==============================================================
8 def loop_solution(tim, prec, values):
9 10001 62502 6.2 0.5 for i, ti in enumerate(tim):
10 10000 1287698 128.8 10.5 x = prec * ti
11 10000 9758633 975.9 79.7 x = np.sin(x)
12 10000 1058995 105.9 8.6 x = np.sum(x)
13 10000 75760 7.6 0.6 values[i] = x
Но помимо фактического времени, затраченного на потребление памяти, есть и другое. Для решения вашей петли требуется память O(n), потому что вы всегда обрабатываете элементы n. Однако для широковещательного решения требуется O(n*n) память. Вероятно, вам придется подождать некоторое время, если вы используете resolution=20000 с вашим циклом, но для него все равно потребуется только 8bytes/element * 20000 element ~= 160kB, но с трансляцией вам понадобится ~3GB. И это пренебрегает постоянным множителем (например, временными массивами или промежуточными массивами)! Предположим, что вы пойдете еще дальше, вы быстро исчерпаете память!
Время для повторного подсчета очков:
- Если вы выполняете цикл python над отдельными элементами в массиве numpy, вы делаете это неправильно.
- Если вы зацикливаете на subarrays из numpy-массива, убедитесь, что служебные данные вызова функции в каждом цикле небрежны по сравнению с временем, проведенным в этой функции.
- Если вы передаете массивы numpy, убедитесь, что у вас не хватает памяти.
Но самая важная точка оптимизации по-прежнему:
-
Только оптимизируйте код, если он слишком медленный! Если он слишком медленный, то оптимизируйте только после профилирования вашего кода.
-
Не слепо доверять упрощенным операциям и никогда не оптимизировать без профилирования.
Последняя мысль:
Такие функции, которые требуют цикла или широковещания, могут быть легко реализованы с использованием numba или numexpr, если их уже нет существующее решение в numpy или scipy.
Например, функция numba, которая сочетает эффективность памяти с решением петли со скоростью широковещательного решения при низком resolutions, будет выглядеть так:
from numba import njit
import math
@njit
def numba_solution(tim, prec, values):
size = tim.size
for i in range(size):
ti = tim[i]
x = 0
for j in range(size):
x += math.sin(prec[j] * ti)
values[i] = x
Как указано в комментариях numexpr, можно также оценить широковещательный расчет очень быстро и без, требующий O(n*n) памяти:
>>> import numexpr
>>> tim_2d = tim[:, np.newaxis]
>>> numexpr.evaluate('sum(sin(tim_2d * prec), axis=1)')