Ответ 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)')