Эффективная оценка функции в каждой ячейке массива NumPy

Учитывая массив A NumPy , каков самый быстрый / наиболее эффективный способ применить одну и ту же функцию f к каждой ячейке?

1. Предположим, что мы присвоим A (i, j) значение f (A (i, j)) .

2. Функция f не имеет двоичного вывода, поэтому операции маски (ing) не помогут.

Является ли "очевидная" итерация двойного цикла (через каждую ячейку) оптимальным решением?

Вы можете просто векторизовать функцию, а затем применять ее непосредственно к массиву Numpy каждый раз, когда она вам нужна:

import numpy as np

def f(x):
    return x * x + 3 * x - 2 if x > 0 else x * 5 + 8

f = np.vectorize(f)  # or use a different name if you want to keep the original f

result_array = f(A)  # if A is your Numpy array

Вероятно, лучше указать явный тип вывода напрямую при векторизации:

f = np.vectorize(f, otypes=[np.float])

Аналогичный вопрос: отображение массива NumPy на месте . Если вы можете найти ufunc для своего f (), вам следует использовать параметр out.

Если вы работаете с числами и f(A(i,j)) = f(A(j,i)), вы можете использовать scipy.spatial.distance.cdist, определяя f как расстояние между A(i)и A(j).

Я считаю, что нашел лучшее решение. Идея изменить функцию на универсальную функцию Python (см. Документацию ), которая может выполнять параллельные вычисления под капотом.

Можно написать свою собственную настройку ufuncна C, что, безусловно, более эффективно, или путем вызова np.frompyfuncвстроенного фабричного метода. После тестирования это более эффективно, чем np.vectorize:

f = lambda x, y: x * y
f_arr = np.frompyfunc(f, 2, 1)
vf = np.vectorize(f)
arr = np.linspace(0, 1, 10000)

%timeit f_arr(arr, arr) # 307ms
%timeit f_arr(arr, arr) # 450ms

Я также протестировал образцы большего размера, и улучшение пропорционально. Для сравнения производительности других методов см. Этот пост

Когда 2d-массив (или nd-массив) является C- или F-смежным, тогда эта задача отображения функции на 2d-массив практически такая же, как задача отображения функции на 1d-массив - мы просто нужно рассматривать это таким образом, например, через np.ravel(A,'K').

Возможное решение для 1d-массива обсуждалось, например, здесь .

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

У Numpy уже есть оборудование для обработки топоров в наилучшем возможном порядке. Одна из возможностей использовать эту технику np.vectorize. Однако в документации numpy np.vectorizeуказано, что она «предоставляется в первую очередь для удобства, а не для производительности» - медленная функция python остается медленной функцией python со всеми связанными накладными расходами! Другая проблема - это огромное потребление памяти - см., Например, этот SO-пост .

Когда кто-то хочет иметь производительность C-функции, но использовать механизм numpy, хорошим решением является использование numba для создания ufunc, например:

# runtime generated C-function as ufunc
import numba as nb
@nb.vectorize(target="cpu")
def nb_vf(x):
    return x+2*x*x+4*x*x*x

Он легко превосходит, np.vectorizeно также, когда та же функция будет выполняться как умножение / сложение массива numpy, т.е.

# numpy-functionality
def f(x):
    return x+2*x*x+4*x*x*x

# python-function as ufunc
import numpy as np
vf=np.vectorize(f)
vf.__name__="vf"

См. Приложение к этому ответу для кода измерения времени:

Версия Numba (зеленая) примерно в 100 раз быстрее, чем функция python (т.е. np.vectorize), что неудивительно. Но это также примерно в 10 раз быстрее, чем numpy-функциональность, потому что версия numbas не требует промежуточных массивов и, таким образом, использует кеш более эффективно.

Хотя подход numba к ufunc - это хороший компромисс между удобством использования и производительностью, это все еще не лучшее, что мы можем сделать. Тем не менее, не существует серебряной пули или наилучшего подхода для любой задачи - нужно понимать, в чем заключаются ограничения и как их можно уменьшить.

Например, для трансцендентных функций (например exp, sin, cos) Numba не дает каких - либо преимуществ по сравнению с NumPy - х np.exp(нет временных массивов , созданных - основной источник ускорения). Однако моя установка Anaconda использует Intel VML для векторов больше 8192 - он просто не может этого сделать, если память не является непрерывной. Поэтому, возможно, лучше скопировать элементы в непрерывную память, чтобы иметь возможность использовать Intel VML:

import numba as nb
@nb.vectorize(target="cpu")
def nb_vexp(x):
    return np.exp(x)

def np_copy_exp(x):
    copy = np.ravel(x, 'K')
    return np.exp(copy).reshape(x.shape) 

Для справедливости сравнения я отключил распараллеливание VML (см. Код в приложении):

Как можно видеть, как только VML начинает работать, накладные расходы на копирование более чем компенсируются. Однако, когда объем данных становится слишком большим для кэша L3, преимущество становится минимальным, поскольку задача снова становится ограниченной пропускной способностью памяти.

С другой стороны, numba также может использовать Intel SVML, как объясняется в этом посте :

from llvmlite import binding
# set before import
binding.set_option('SVML', '-vector-library=SVML')

import numba as nb

@nb.vectorize(target="cpu")
def nb_vexp_svml(x):
    return np.exp(x)

а использование VML с распараллеливанием дает:

Версия numba имеет меньше накладных расходов, но для некоторых размеров VML превосходит SVML, даже несмотря на дополнительные накладные расходы на копирование, что неудивительно, поскольку ufuncs numba не распараллеливаются.

Тэг:

A. сравнение полиномиальной функции:

import perfplot
perfplot.show(
    setup=lambda n: np.random.rand(n,n)[::2,::2],
    n_range=[2**k for k in range(0,12)],
    kernels=[
        f,
        vf, 
        nb_vf
        ],
    logx=True,
    logy=True,
    xlabel='len(x)'
    ) 

Б. сравнение exp:

import perfplot
import numexpr as ne # using ne is the easiest way to set vml_num_threads
ne.set_vml_num_threads(1)
perfplot.show(
    setup=lambda n: np.random.rand(n,n)[::2,::2],
    n_range=[2**k for k in range(0,12)],
    kernels=[
        nb_vexp, 
        np.exp,
        np_copy_exp,
        ],
    logx=True,
    logy=True,
    xlabel='len(x)',
    )

Все приведенные выше ответы хорошо сравниваются, но если вам нужно использовать настраиваемую функцию для сопоставления, а у вас есть numpy.ndarray, и вам нужно сохранить форму массива.

У меня сравните всего два, но форма сохранит ndarray. Для сравнения я использовал массив с 1 миллионом записей. Здесь я использую квадратную функцию. Я представляю общий случай для n-мерного массива. Для двухмерного просто сделайте iterдля 2D.

import numpy, time

def A(e):
    return e * e

def timeit():
    y = numpy.arange(1000000)
    now = time.time()
    numpy.array([A(x) for x in y.reshape(-1)]).reshape(y.shape)        
    print(time.time() - now)
    now = time.time()
    numpy.fromiter((A(x) for x in y.reshape(-1)), y.dtype).reshape(y.shape)
    print(time.time() - now)
    now = time.time()
    numpy.square(y)  
    print(time.time() - now)

Вывод

>>> timeit()
1.162431240081787    # list comprehension and then building numpy array
1.0775556564331055   # from numpy.fromiter
0.002948284149169922 # using inbuilt function

здесь вы можете четко видеть numpy.fromiterфункцию квадрата пользователя, используйте любую на ваш выбор. Если ваша функция зависит от i, j индексов массива, повторяйте размер массива, например for ind in range(arr.size), используйте, numpy.unravel_indexчтобы получить i, j, ..на основе вашего 1D-индекса и формы массива numpy.unravel_index

Эти ответы вдохновлены моим ответом на другой вопрос здесь