Comment vectoriser série de Fourier somme partielle en numpy

Question

Compte tenu des coefficients de séries de Fourier a[n] et b[n] (par cosinus et sinus respectivement) d'une fonction par période T et t un intervalle également espacés le code suivant permettra d'évaluer la somme partielle pour tous les points dans l'intervalle t (a, b, t sont tous des tableaux numpy). Il est précisé que len (t) <> len (a). Ma question est: Est-ce que cette boucle for peut être vectorisée?

Answer 1

est ici une approche vectorisé utilisant broadcasting pour créer la version tableau 2D de cosinus/entrée sinus: 2*pi*n*t/T puis en utilisant matrix-multiplication avec np.dot pour le sum-reduction -

r = np.arange(1,len(a)) 
S = 2*np.pi*r[:,None]*t/T 
cS = np.cos(S) 
sS = np.sin(S) 
out = a[1:].dot(cS) - b[1:].dot(sS) + a[0] 

plus poussée de la performance

Pour plus de boost, nous pouvons utiliser numexpr module pour calculer ces étapes trignométriques -

import numexpr as ne 
cS = ne.evaluate('cos(S)') 
sS = ne.evaluate('sin(S)') 

Test d'exécution -

approches -

def original_app(t,a,b,T): 
    yn=np.ones(len(t))*a[0] 
    for n in range(1,len(a)): 
     yn=yn+(a[n]*np.cos(2*np.pi*n*t/T)-b[n]*np.sin(2*np.pi*n*t/T)) 
    return yn 

def vectorized_app(t,a,b,T):  
    r = np.arange(1,len(a)) 
    S = (2*np.pi/T)*r[:,None]*t 
    cS = np.cos(S) 
    sS = np.sin(S) 
    return a[1:].dot(cS) - b[1:].dot(sS) + a[0] 

def vectorized_app_v2(t,a,b,T):  
    r = np.arange(1,len(a)) 
    S = (2*np.pi/T)*r[:,None]*t 
    cS = ne.evaluate('cos(S)') 
    sS = ne.evaluate('sin(S)') 
    return a[1:].dot(cS) - b[1:].dot(sS) + a[0] 

également, y compris la fonction PP de poste de Panzer @ Paul.

synchronisations -

In [22]: # Setup inputs 
    ...: n = 10000 
    ...: t = np.random.randint(0,9,(n)) 
    ...: a = np.random.randint(0,9,(n)) 
    ...: b = np.random.randint(0,9,(n)) 
    ...: T = 3.45 
    ...: 

In [23]: print np.allclose(original_app(t,a,b,T), vectorized_app(t,a,b,T)) 
    ...: print np.allclose(original_app(t,a,b,T), vectorized_app_v2(t,a,b,T)) 
    ...: print np.allclose(original_app(t,a,b,T), PP(t,a,b,T)) 
    ...: 
True 
True 
True 

In [25]: %timeit original_app(t,a,b,T) 
    ...: %timeit vectorized_app(t,a,b,T) 
    ...: %timeit vectorized_app_v2(t,a,b,T) 
    ...: %timeit PP(t,a,b,T) 
    ...: 
1 loops, best of 3: 6.49 s per loop 
1 loops, best of 3: 6.24 s per loop 
1 loops, best of 3: 1.54 s per loop 
1 loops, best of 3: 1.96 s per loop 

Answer 2

ne peut pas battre numexpr, mais si ce n'est pas disponible, nous pouvons économiser sur les transcendantaux (tests et le code étalonnage fortement basé sur @ code de Divakar au cas où vous ne l'avez pas remarqué; -)):

import numpy as np 
from timeit import timeit 

def PP(t,a,b,T): 
    CS = np.empty((len(t), len(a)-1), np.complex) 
    CS[...] = np.exp(2j*np.pi*(t[:, None])/T) 
    np.cumprod(CS, axis=-1, out=CS) 
    return a[1:].dot(CS.T.real) - b[1:].dot(CS.T.imag) + a[0] 

def original_app(t,a,b,T): 
    yn=np.ones(len(t))*a[0] 
    for n in range(1,len(a)): 
     yn=yn+(a[n]*np.cos(2*np.pi*n*t/T)-b[n]*np.sin(2*np.pi*n*t/T)) 
    return yn 

def vectorized_app(t,a,b,T):  
    r = np.arange(1,len(a)) 
    S = 2*np.pi*r[:,None]*t/T 
    cS = np.cos(S) 
    sS = np.sin(S) 
    return a[1:].dot(cS) - b[1:].dot(sS) + a[0] 

n = 1000 
t = 2000 
t = np.random.randint(0,9,(t)) 
a = np.random.randint(0,9,(n)) 
b = np.random.randint(0,9,(n)) 
T = 3.45 


print(np.allclose(original_app(t,a,b,T), vectorized_app(t,a,b,T))) 
print(np.allclose(original_app(t,a,b,T), PP(t,a,b,T))) 

print('{:18s} {:9.6f}'.format('orig', timeit(lambda: original_app(t,a,b,T), number=10)/10)) 
print('{:18s} {:9.6f}'.format('Divakar no numexpr', timeit(lambda: vectorized_app(t,a,b,T), number=10)/10)) 
print('{:18s} {:9.6f}'.format('PP', timeit(lambda: PP(t,a,b,T), number=10)/10)) 

Prints:

True 
True 
orig    0.166903 
Divakar no numexpr 0.179617 
PP     0.060817 

btw. si delta t divise T on peut potentiellement économiser plus, ou même exécuter le plein fft et jeter ce qui est trop.

Answer 3

Ce n'est pas vraiment une autre réponse mais un commentaire sur celui de @Paul Panzer, écrit en réponse car j'avais besoin de poster du code. S'il y a un moyen d'afficher du code correctement formaté dans un commentaire, veuillez le conseiller.

Inspiré par Panzer idée cumprod @ Paul, je suis venu avec les éléments suivants:

an = ones((len(a)-1,len(te)))*2j*pi*te/T 
CS = exp(cumsum(an,axis=0)) 
out = (a[1:].dot(CS.real) - b[1:].dot(CS.imag)) + a[0] 

Bien qu'il semble bien vectorisé et produit des résultats corrects, sa performance est misérable. Il est non seulement beaucoup plus lent que le cumprod, ce qui est attendu comme len(a)-1 exponentiations plus sont faites, mais 50% plus lent que la version originale nonvectorisée. Quelle est la cause de cette mauvaise performance?