Можно ли векторизовать рекурсивный расчет массива NumPy, где каждый элемент зависит от предыдущего?
T(i) = Tm(i) + (T(i-1)-Tm(i))**(-tau(i))
Tm а также tau являются векторами NumPy той же длины, которые были рассчитаны ранее, и желательно создать новый вектор T, i включен только для указания индекса элемента для того, что требуется.
Необходим ли цикл for для этого случая?
5 ответов
Вы можете подумать, что это будет работать:
import numpy as np
n = len(Tm)
t = np.empty(n)
t[0] = 0 # or whatever the initial condition is
t[1:] = Tm[1:] + (t[0:n-1] - Tm[1:])**(-tau[1:])
но это не так: на самом деле вы не можете сделать рекурсию в numpy таким образом (поскольку numpy вычисляет всю RHS и затем назначает ее LHS).
Поэтому, если вы не можете придумать нерекурсивную версию этой формулы, вы застряли с явным циклом:
tt = np.empty(n)
tt[0] = 0.
for i in range(1,n):
tt[i] = Tm[i] + (tt[i-1] - Tm[i])**(-tau[i])
Я выполнил некоторые тесты, и в 2018 году использование Numba - это первый вариант, который люди должны пытаться ускорить рекурсивными функциями в Numpy (предложение Aronstef). Numba уже предустановлена в пакете Anaconda и имеет один из самых быстрых раз (примерно в 20 раз быстрее, чем любой Python). В 2018 году Python поддерживает аннотации @numba без дополнительных действий (как минимум версии 3.6 и 3.7). Вот два контрольных показателя: один выполнен на 2018-10-20, а другой на 2016-05-18.
И, как упоминал Джаффе, в 2018 году все еще невозможно векторизовать рекурсивные функции. Я проверил векторизацию Aronstef, и она не работает.
Бенчмарки отсортированы по времени выполнения:
-----------------------------------
|Variant |2018-10 |2016-05 |
-----------------------------------
|Pure C | na | 2.75 ms|
|C extension | na | 6.22 ms|
|Cython float32 | 1.01 ms| na |
|Cython float64 | 1.05 ms| 6.26 ms|
|Fortran f2py | na | 6.78 ms|
|Numba float32 | 2.81 ms| na |
|(Aronstef) | | |
|Numba float64 | 5.28 ms| na |
|Append to list |48.2 ms|91.0 ms|
|Using a.item() |58.3 ms|74.4 ms|
|np.fromiter() |60.0 ms|78.1 ms|
|Loop over Numpy|71.9 ms|87.9 ms|
|(Jaffe) | | |
|Loop over Numpy|74.4 ms| na |
|(Aronstef) | | |
-----------------------------------
Соответствующий код указан в конце ответа.
Я не проверял Pure C в 2018 году, но полагаю, что он все еще самый быстрый на основе предыдущего теста.
Я также не проверял расширение C в 2018 году и думаю, что оно почти такое же время, как и Cython на основе предыдущего теста.
Fortran было очень сложно отлаживать и компилировать, поэтому я не проверял версию f2py в 2018 году. И все равно она была хуже, чем Cython.
У меня есть следующие настройки в 2018 году:
Processor: Intel i7-7500U 2.7GHz
Versions:
Python: 3.7.0
Numba: 0.39.0
Cython: 0.28.5
Numpy: 1.15.1
Рекомендуемый код Numba с использованием float32 (от Aronstef):
@numba.jit("float32[:](float32[:], float32[:])", nopython=False, nogil=True)
def calc_py_jit32(Tm_, tau_):
tt = np.empty(len(Tm_),dtype="float32")
tt[0] = Tm_[0]
for i in range(1, len(Tm_)):
tt[i] = Tm_[i] - (tt[i-1] + Tm_[i])**(-tau_[i])
return tt[1:]
Весь другой код:
Создание данных (например, комментарий Aronstef + Mike T):
np.random.seed(0)
n = 100000
Tm = np.cumsum(np.random.uniform(0.1, 1, size=n).astype('float64'))
tau = np.random.uniform(-1, 0, size=n).astype('float64')
ar = np.column_stack([Tm,tau])
Tm32 = Tm.astype('float32')
tau32 = tau.astype('float32')
Tm_l = list(Tm)
tau_l = list(tau)
Код в 2016 году немного отличался, поскольку я использовал функцию abs() для предотвращения nans, а не вариант Mike T. В 2018 году функция точно такая же, как написал OP (Original Poster).
Cython float32, использующий магию Jupyter %%. Функция может быть использована непосредственно в Python, Cython нужен компилятор C++, в котором был скомпилирован Python. Установка правильной версии компилятора Visual C++ (для Windows) может быть проблематичной:
%%cython
import cython
import numpy as np
cimport numpy as np
from numpy cimport ndarray
cdef extern from "math.h":
np.float32_t exp(np.float32_t m)
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.infer_types(True)
@cython.initializedcheck(False)
def cy_loop32(np.float32_t[:] Tm,np.float32_t[:] tau,int alen):
cdef np.float32_t[:] T=np.empty(alen, dtype=np.float32)
cdef int i
T[0]=0.0
for i in range(1,alen):
T[i] = Tm[i] + (T[i-1] - Tm[i])**(-tau[i])
return T
Cython float64 с использованием Jupyter %% magic. Функция может быть использована непосредственно в Python:
%%cython
cdef extern from "math.h":
double exp(double m)
import cython
import numpy as np
cimport numpy as np
from numpy cimport ndarray
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.infer_types(True)
@cython.initializedcheck(False)
def cy_loop(double[:] Tm,double[:] tau,int alen):
cdef double[:] T=np.empty(alen)
cdef int i
T[0]=0.0
for i in range(1,alen):
T[i] = Tm[i] + (T[i-1] - Tm[i])**(-tau[i])
return T
Numba float64:
@numba.jit("float64[:](float64[:], float64[:])", nopython=False, nogil=True)
def calc_py_jit(Tm_, tau_):
tt = np.empty(len(Tm_),dtype="float64")
tt[0] = Tm_[0]
for i in range(1, len(Tm_)):
tt[i] = Tm_[i] - (tt[i-1] + Tm_[i])**(-tau_[i])
return tt[1:]
Добавить в список Самое быстрое не скомпилированное решение:
def rec_py_loop(Tm,tau,alen):
T = [Tm[0]]
for i in range(1,alen):
T.append(Tm[i] - (T[i-1] + Tm[i])**(-tau[i]))
return np.array(T)
Использование a.item ():
def rec_numpy_loop_item(Tm_,tau_):
n_ = len(Tm_)
tt=np.empty(n_)
Ti=tt.item
Tis=tt.itemset
Tmi=Tm_.item
taui=tau_.item
Tis(0,Tm_[0])
for i in range(1,n_):
Tis(i,Tmi(i) - (Ti(i-1) + Tmi(i))**(-taui(i)))
return tt[1:]
np.fromiter ():
def it(Tm,tau):
T=Tm[0]
i=0
while True:
yield T
i+=1
T=Tm[i] - (T + Tm[i])**(-tau[i])
def rec_numpy_iter(Tm,tau,alen):
return np.fromiter(it(Tm,tau), np.float64, alen)[1:]
Цикл по Numpy (основанный на идее Jaffe):
def rec_numpy_loop(Tm,tau,alen):
tt=np.empty(alen)
tt[0]=Tm[0]
for i in range(1,alen):
tt[i] = Tm[i] - (tt[i-1] + Tm[i])**(-tau[i])
return tt[1:]
Зацикливайтесь на Numpy (код Аронстефа). На моем компьютере float64 тип по умолчанию для np.empty,
def calc_py(Tm_, tau_):
tt = np.empty(len(Tm_),dtype="float64")
tt[0] = Tm_[0]
for i in range(1, len(Tm_)):
tt[i] = (Tm_[i] - (tt[i-1] + Tm_[i])**(-tau_[i]))
return tt[1:]
Чистый C без использования Python совсем. Версия от 2016 года (с функцией fabs()):
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <windows.h>
#include <sys\timeb.h>
double randn() {
double u = rand();
if (u > 0.5) {
return sqrt(-1.57079632679*log(1.0 - pow(2.0 * u - 1, 2)));
}
else {
return -sqrt(-1.57079632679*log(1.0 - pow(1 - 2.0 * u,2)));
}
}
void rec_pure_c(double *Tm, double *tau, int alen, double *T)
{
for (int i = 1; i < alen; i++)
{
T[i] = Tm[i] + pow(fabs(T[i - 1] - Tm[i]), (-tau[i]));
}
}
int main() {
int N = 100000;
double *Tm= calloc(N, sizeof *Tm);
double *tau = calloc(N, sizeof *tau);
double *T = calloc(N, sizeof *T);
double time = 0;
double sumtime = 0;
for (int i = 0; i < N; i++)
{
Tm[i] = randn();
tau[i] = randn();
}
LARGE_INTEGER StartingTime, EndingTime, ElapsedMicroseconds;
LARGE_INTEGER Frequency;
for (int j = 0; j < 1000; j++)
{
for (int i = 0; i < 3; i++)
{
QueryPerformanceFrequency(&Frequency);
QueryPerformanceCounter(&StartingTime);
rec_pure_c(Tm, tau, N, T);
QueryPerformanceCounter(&EndingTime);
ElapsedMicroseconds.QuadPart = EndingTime.QuadPart - StartingTime.QuadPart;
ElapsedMicroseconds.QuadPart *= 1000000;
ElapsedMicroseconds.QuadPart /= Frequency.QuadPart;
if (i == 0)
time = (double)ElapsedMicroseconds.QuadPart / 1000;
else {
if (time > (double)ElapsedMicroseconds.QuadPart / 1000)
time = (double)ElapsedMicroseconds.QuadPart / 1000;
}
}
sumtime += time;
}
printf("1000 loops,best of 3: %.3f ms per loop\n",sumtime/1000);
free(Tm);
free(tau);
free(T);
}
Фортран ф2пи. Функция может быть использована с Python, Версия от 2016 года (с функцией abs()):
subroutine rec_fortran(tm,tau,alen,result)
integer*8, intent(in) :: alen
real*8, dimension(alen), intent(in) :: tm
real*8, dimension(alen), intent(in) :: tau
real*8, dimension(alen) :: res
real*8, dimension(alen), intent(out) :: result
res(1)=0
do i=2,alen
res(i) = tm(i) + (abs(res(i-1) - tm(i)))**(-tau(i))
end do
result=res
end subroutine rec_fortran
В некоторых случаях возможна такая рекурсия, а именно, если для формулы рекурсии есть Ufunc NumPy, например
c = numpy.arange(10.)
numpy.add(c[:-1], c[1:], c[1:])
Это вычисляет накопленные суммы за c на месте, используя выходной параметр numpy.add, Это не может быть написано как
c[1:] = c[:-1] + c[1:]
потому что теперь результат добавления является временным, который копируется в c[1:] после того, как вычисление закончено.
Самое естественное, что можно попробовать сейчас, это определить свой собственный ufunc:
def f(T, Tm, tau):
return Tm + (T - Tm)**(-tau)
uf = numpy.frompyfunc(f, 3, 1)
Но по причинам, которые мне недоступны, следующее не работает:
uf(T[:-1], Tm[1:], tau[1:], T[1:])
Видимо, результат не написан напрямую T[1:], но хранится во временном хранилище и копируется после завершения операции. Даже если бы это сработало, я бы не ожидал от этого ускорения по сравнению с обычным циклом, поскольку для каждой записи необходимо вызывать функцию Python.
Если вы действительно хотите избежать цикла Python, вам, вероятно, нужно перейти на Cython или ctypes.
Обновление: 21-10-2018 Я исправил свой ответ на основе комментариев.
Можно векторизовать операции над векторами, если расчет не является рекурсивным. Поскольку рекурсивная операция зависит от предыдущего вычисленного значения, параллельная обработка операции невозможна. Поэтому это не работает:
def calc_vect(Tm_, tau_):
return Tm_[1:] - (Tm_[:-1] + Tm_[1:]) ** (-tau_[1:])
Поскольку (последовательная обработка / цикл) необходима, наилучшая производительность достигается при максимально возможном приближении к оптимизированному машинному коду, поэтому Numba и Cython - лучшие ответы здесь.
Подход Numba может быть достигнут следующим образом:
init_string = """
from math import pow
import numpy as np
from numba import jit, float32
np.random.seed(0)
n = 100000
Tm = np.cumsum(np.random.uniform(0.1, 1, size=n).astype('float32'))
tau = np.random.uniform(-1, 0, size=n).astype('float32')
def calc_python(Tm_, tau_):
tt = np.empty(len(Tm_))
tt[0] = Tm_[0]
for i in range(1, len(Tm_)):
tt[i] = Tm_[i] - pow(tt[i-1] + Tm_[i], -tau_[i])
return tt
@jit(float32[:](float32[:], float32[:]), nopython=False, nogil=True)
def calc_numba(Tm_, tau_):
tt = np.empty(len(Tm_))
tt[0] = Tm_[0]
for i in range(1, len(Tm_)):
tt[i] = Tm_[i] - pow(tt[i-1] + Tm_[i], -tau_[i])
return tt
"""
import timeit
py_time = timeit.timeit('calc_python(Tm, tau)', init_string, number=100)
numba_time = timeit.timeit('calc_numba(Tm, tau)', init_string, number=100)
print("Python Solution: {}".format(py_time))
print("Numba Soltution: {}".format(numba_time))
Сравнение времени выполнения функций Python и Numba:
Python Solution: 54.58057559299999
Numba Soltution: 1.1389029540000024
Это хороший вопрос. Мне также интересно узнать, возможно ли это, но пока я не нашел способа сделать это, за исключением некоторых простых случаев.
Опция 1. numpy.ufunc.accumulate
Как отметил @Karl Knechtel, этот вариант кажется многообещающим. Вам нужно создать ufunc первый. Эта веб-страница объясняет, как это сделать.
В простом случае рекуррентной функции, которая принимает два скаляра в качестве входных данных и выводит один масштабатор, кажется, что это работает:
import numpy as np
def test_add(x, data):
return x + data
assert test_add(1, 2) == 3
assert test_add(2, 3) == 5
# Make a Numpy ufunc from my test_add function
test_add_ufunc = np.frompyfunc(t_next, 2, 1)
assert test_add_ufunc(1, 2) == 3
assert test_add_ufunc(2, 3) == 5
assert np.all(test_add_ufunc([1, 2], [2, 3]) == [3, 5])
data_sequence = np.array([1, 2, 3, 4])
f_out = test_add_ufunc.accumulate(data_sequence, dtype=object)
assert np.array_equal(f_out, [1, 3, 6, 10])
[Обратите внимание dtype=object аргумент, который необходим, как описано на веб-странице, указанной выше].
Но в вашем случае (и моем) мы хотим вычислить рекуррентное уравнение, которое имеет более одного ввода данных (и, возможно, более одной переменной состояния).
Когда я попробовал это с помощью ufunc.accumulate подход выше я получил ValueError: accumulate only supported for binary functions.
Если кто-нибудь знает способ обойти это ограничение, мне было бы очень интересно.
Вариант 2. Встроенная функция накопления в Python
Между тем, это решение не совсем соответствует тому, что вы хотели с точки зрения векторизованного вычисления в numpy, но оно, по крайней мере, позволяет избежать цикла for.
from itertools import accumulate, chain
def t_next(t, data):
Tm, tau = data # Unpack more than one data input
return Tm + (t - Tm)**tau
assert t_next(2, (0.38, 0)) == 1.38
t0 = 2 # Initial t
Tm_values = np.array([0.38, 0.88, 0.56, 0.67, 0.45, 0.98, 0.58, 0.72, 0.92, 0.82])
tau_values = np.linspace(0, 0.9, 10)
# Combine the input data into a 2D array
data_sequence = np.vstack([Tm_values, tau_values]).T
t_out = np.fromiter(accumulate(chain([t0], data_sequence), t_next), dtype=float)
print(t_out)
# [2. 1.38 1.81303299 1.60614649 1.65039964 1.52579703
# 1.71878078 1.66109554 1.67839293 1.72152195 1.73091672]
# Slightly more readable version possible in Python 3.8+
t_out = np.fromiter(accumulate(data_sequence, t_next, initial=t0), dtype=float)
print(t_out)
# [2. 1.38 1.81303299 1.60614649 1.65039964 1.52579703
# 1.71878078 1.66109554 1.67839293 1.72152195 1.73091672]
Чтобы основываться на ответе NPE, я согласен, что где-то должен быть цикл. Возможно, ваша цель - избежать накладных расходов, связанных с циклом Python for? В этом случае numpy.fromiter обходит цикл for, но только немного:
Используя очень простое отношение рекурсии,
x[i+1] = x[i] + 0.1
я получил
#FOR LOOP
def loopit(n):
x = [0.0]
for i in range(n-1): x.append(x[-1] + 0.1)
return np.array(x)
#FROMITER
#define an iterator (a better way probably exists -- I'm a novice)
def it():
x = 0.0
while True:
yield x
x += 0.1
#use the iterator with np.fromiter
def fi_it(n):
return np.fromiter(it(), np.float, n)
%timeit -n 100 loopit(100000)
#100 loops, best of 3: 31.7 ms per loop
%timeit -n 100 fi_it(100000)
#100 loops, best of 3: 18.6 ms per loop
Интересно, что предварительное выделение пустого массива приводит к существенной потере производительности. Для меня это загадка, хотя я предполагаю, что при доступе к элементу массива должно быть больше накладных расходов, чем при добавлении в список.
def loopit(n):
x = np.zeros(n)
for i in range(n-1): x[i+1] = x[i] + 0.1
return x
%timeit -n 100 loopit(100000)
#100 loops, best of 3: 50.1 ms per loop