Можно ли векторизовать рекурсивное вычисление массива NumPy, где каждый элемент зависит от предыдущего?

Вы могли подумать, что это сработает:

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])

Обновление 2019. Код Numba был нарушен с новой версией Numba. Замена dtype="float32" на dtype=np.float32 решила эту проблему.

Я выполнил несколько тестов, и в 2019 году использование Numba - это первый вариант, который люди должны попытаться ускорить рекурсивные функции в Numpy (скорректированное предложение Аронстефа). Numba уже предустановлена ​​в пакете Anaconda и имеет одно из самых быстрых времен (примерно в 20 раз быстрее, чем любой Python). В 2019 Python поддерживает аннотации @numba без дополнительных действий (как минимум версии 3.6, 3.7 и 3.8). Вот три теста: выполнено 5 декабря 2019 г., 20 октября 2018 г. и 18 мая 2016 г.

И, как упоминал Яффе, в 2018 году все еще невозможно векторизовать рекурсивные функции. Я проверил векторизацию Aronstef, и она НЕ работает.

Бенчмарки, отсортированные по времени выполнения:

-------------------------------------------
|Variant        |2019-12 |2018-10 |2016-05 |
-------------------------------------------
|Pure C         |   na   |   na   | 2.75 ms|
|C extension    |   na   |   na   | 6.22 ms|
|Cython float32 | 0.55 ms| 1.01 ms|   na   |
|Cython float64 | 0.54 ms| 1.05 ms| 6.26 ms|
|Fortran f2py   | 4.65 ms|   na   | 6.78 ms|
|Numba float32  |73.0  ms| 2.81 ms|   na   |
|(Aronstef)     |        |        |        |
|Numba float32v2| 1.82 ms| 2.81 ms|   na   |
|Numba float64  |78.9  ms| 5.28 ms|   na   |
|Numba float64v2| 4.49 ms| 5.28 ms|   na   |
|Append to list |73.3  ms|48.2  ms|91.0  ms|
|Using a.item() |36.9  ms|58.3  ms|74.4  ms|
|np.fromiter()  |60.8  ms|60.0  ms|78.1  ms|
|Loop over Numpy|71.3  ms|71.9  ms|87.9  ms|
|(Jaffe)        |        |        |        |
|Loop over Numpy|74.6  ms|74.4  ms|   na   |
|(Aronstef)     |        |        |        |
-------------------------------------------

Соответствующий код указан в конце ответа.

Кажется, что со временем времена Numba и Cython поправляются. Теперь оба они быстрее, чем Fortran f2py. Cython теперь быстрее в 8,6 раза, а Numba 32bit быстрее в 2,5 раза. В 2016 году Фортран было очень сложно отлаживать и компилировать. Так что теперь нет причин использовать Фортран вообще.

Я не проверял расширения Pure C и C в 2019 и 2018 годах, потому что их непросто скомпилировать в ноутбуках Jupyter.

В 2019 году у меня была следующая установка:

Processor: Intel i5-9600K 3.70GHz
Versions:
Python:  3.8.0
Numba:  0.46.0
Cython: 0.29.14
Numpy:  1.17.4

В 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 (скорректированный Аронстеф):

@numba.jit("float32[:](float32[:], float32[:])", nopython=True, nogil=True)
def calc_py_jit32v2(Tm_, tau_):
    tt = np.empty(len(Tm_),dtype=np.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, а не вариант Майка Т. В 2018 году функция точно такая же, как написал OP (Original Poster).

Cython float32 с использованием Jupyter %% magic. Функцию можно использовать непосредственно в 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_jitv2(Tm_, tau_):
    tt = np.empty(len(Tm_),dtype=np.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 (на основе идеи Джеффа):

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);
}

Fortran f2py. Функцию можно использовать с 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

Обновление: 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))

Сравнение Timeit функций Python и Numba:

Python Solution: 54.58057559299999
Numba Soltution: 1.1389029540000024

Чтобы основываться на ответе NPE, я согласен, что где-то должен быть цикл. Возможно, ваша цель - избежать накладных расходов, связанных с циклом for Python? В этом случае 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

Интересно, что предварительное выделение массива numpy приводит к существенной потере производительности. Для меня это загадка, хотя я предполагаю, что с доступом к элементу массива должно быть больше накладных расходов, чем с добавлением в список.

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

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

Вариант 1. numpy.ufunc.accumulate

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(test_add, 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])

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]