Est-il possible de vectoriser le calcul récursif d'un tableau NumPy où chaque élément dépend du précédent ?

Question

T(i) = Tm(i) + (T(i-1)-Tm(i))**(-tau(i))

Tm y tau sont des vecteurs NumPy de même longueur qui ont été calculés précédemment, et le désir est de créer un nouveau vecteur T . Le site i est inclus uniquement pour indiquer l'index de l'élément pour ce qui est souhaité.

Une boucle for est-elle nécessaire dans ce cas ?

Answer 1

On pourrait penser que ça marche :

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

mais ce n'est pas le cas : vous ne pouvez pas réellement faire de récursion dans numpy de cette façon (puisque numpy calcule la totalité de la RHS et l'assigne ensuite à la LHS).

Donc, à moins que vous ne trouviez une version non récursive de cette formule, vous êtes coincé avec une boucle explicite :

tt = np.empty(n)
tt[0] = 0.
for i in range(1,n):
    tt[i] = Tm[i] + (tt[i-1] - Tm[i])**(-tau[i])

Answer 2

Mise à jour de 2019. Le code Numba s'est cassé avec la nouvelle version de numba. Changer dtype="float32" a dtype=np.float32 l'a résolu.

J'ai effectué quelques tests et en 2019, en utilisant Numba est la première option que les gens devraient essayer pour accélérer les fonctions récursives dans Numpy (proposition ajustée d'Aronstef). Numba est déjà préinstallé dans le paquetage Anaconda et a l'un des temps les plus rapides (environ 20 fois plus rapide que n'importe quel Python). En 2019, Python prend en charge les annotations @numba sans étapes supplémentaires (au moins les versions 3.6, 3.7 et 3.8). Voici trois benchmarks : réalisés les 2019-12-05, 2018-10-20 et 2016-05-18.

Et, comme mentionné par Jaffe, en 2018 il n'est toujours pas possible de vectoriser les fonctions récursives. J'ai vérifié la vectorisation par Aronstef et cela ne fonctionne PAS.

Benchmarks triés par temps d'exécution :

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

Le code correspondant est fourni à la fin de la réponse.

Il semble qu'avec le temps Numba et Cython s'améliorent. Maintenant, les deux sont plus rapides que Fortran f2py. Cython est 8,6 fois plus rapide maintenant et Numba 32bit est 2,5 fois plus rapide. Fortran était très difficile à déboguer et à compiler en 2016. Donc maintenant, il n'y a aucune raison d'utiliser Fortran du tout.

Je n'ai pas vérifié Pure C et C extension en 2019 et 2018, car il n'est pas facile de les compiler dans les carnets Jupyter.

J'avais la configuration suivante en 2019 :

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

J'avais la configuration suivante en 2018 :

Processor: Intel i7-7500U 2.7GHz
Versions:
Python:  3.7.0
Numba:  0.39.0
Cython: 0.28.5
Numpy:  1.15.1

La recommandation Numba code utilisant float32 (Aronstef ajusté) :

@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:]

Tous les autres codes :

Création de données (comme Aronstef + commentaire de 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)

Le code en 2016 était légèrement différent car j'utilisais la fonction abs() pour prévenir les nans et non la variante de Mike T. En 2018, la fonction est exactement la même que celle écrite par OP (Original Poster).

Cython float32 en utilisant la magie de Jupyter %%. La fonction peut être utilisée directement dans Python . Cython a besoin d'un compilateur C++ dans lequel Python a été compilé. L'installation de la bonne version du compilateur Visual C++ (pour Windows) peut être problématique :

%%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 en utilisant la magie de Jupyter %%. La fonction peut être utilisée directement dans 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:]

Ajouter à la liste . Solution non compilée la plus rapide :

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)

Utilisation de 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:]

Boucle sur Numpy (basé sur l'idée de 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:]

Boucle sur Numpy (code d'Aronstef). Sur mon ordinateur float64 est le type par défaut pour 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:]

Pure C sans utiliser Python du tout. Version de l'année 2016 (avec la fonction 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. La fonction peut être utilisée à partir de Python . Version de l'année 2016 (avec la fonction 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

Answer 3

Mise à jour : 21-10-2018 J'ai corrigé ma réponse en fonction des commentaires.

Il est possible de vectoriser des opérations sur des vecteurs tant que le calcul n'est pas récursif. Comme une opération récursive dépend de la valeur calculée précédemment, il n'est pas possible de traiter l'opération en parallèle. Cela ne fonctionne donc pas :

def calc_vect(Tm_, tau_):
    return Tm_[1:] - (Tm_[:-1] + Tm_[1:]) ** (-tau_[1:])

Puisque (traitement en série / une boucle) est nécessaire, les meilleures performances sont obtenues en se rapprochant le plus possible du code machine optimisé, donc Numba et Cython sont les meilleures réponses ici.

Une approche Numba peut être réalisée comme suit :

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

Comparaison en temps réel des fonctions Python et Numba :

Python Solution: 54.58057559299999
Numba Soltution: 1.1389029540000024

Answer 4

C'est une bonne question. J'aimerais également savoir si cela est possible, mais jusqu'à présent, je n'ai pas trouvé de moyen de le faire, sauf dans certains cas simples.

Option 1. numpy.ufunc.accumulate

Cela semble être une option prometteuse comme mentionné par @Karl Knechtel. Vous devez créer un ufunc d'abord. Cette page web explique comment.

Dans le cas simple d'une fonction récurrente qui prend deux scalaires en entrée et sort un scalaire, cela semble fonctionner :

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

[Notez le dtype=object argument qui est nécessaire, comme expliqué sur la page web liée ci-dessus].

Mais dans votre cas (et le mien), nous voulons calculer une équation récurrente qui a plus d'une entrée de données (et potentiellement plus d'une variable d'état aussi).

Lorsque j'ai essayé de le faire en utilisant le ufunc.accumulate l'approche ci-dessus j'ai obtenu ValueError: accumulate only supported for binary functions .

Si quelqu'un connaît un moyen de contourner cette contrainte, je serais très intéressé.

Option 2. La fonction intégrée de Python accumuler fonction

En attendant, cette solution n'atteint pas tout à fait ce que vous souhaitiez en termes de calcul vectoriel dans numpy, mais elle permet au moins d'éviter une boucle 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]

Answer 5

Pour compléter la réponse de NPE, je suis d'accord pour dire qu'il doit y avoir une boucle quelque part. Peut-être que votre objectif est d'éviter la surcharge associée à une boucle for Python ? Dans ce cas, numpy.fromiter est plus efficace qu'une boucle for, mais seulement de peu :

En utilisant la relation de récursion très simple,

x[i+1] = x[i] + 0.1

Je reçois

#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

Il est intéressant de noter que la pré-affectation d'un tableau numpy entraîne une perte substantielle des performances. C'est un mystère pour moi, bien que je suppose qu'il doit y avoir plus de surcharge associée à l'accès à un élément de tableau qu'à l'ajout à une liste.

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