如何向量化傅里葉級數部分和在numpy的

Question

給定一個功能的具有周期T和t的等距間隔的傅立葉級數的係數（分別爲餘弦和正弦）a[n]和b[n]下面的代碼將評估所有點的部分和在區間t（a,b,t都是numpy陣列）。有人澄清，len（t）<> len（a）。

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

我的問題是：這個for循環是否可以被矢量化？

Answer 1

這裏有一個量化的方法利用broadcasting創建餘弦/正弦輸入的2D陣列版本：np.dot2*pi*n*t/T，然後使用matrix-multiplication爲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] 

進一步的性能提升

爲了進一步提升，我們可以利用numexpr module來計算那些三角函數步 -

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

運行測試 -

途徑 -

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] 

此外，包括@保羅裝甲的帖子功能PP。

計時 -

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

不能擊敗numexpr，但如果它不提供，我們可以節省上超越數（很大程度上基於@測試和基準測試代碼Divakar在你沒有注意到的情況下代碼; - ））：

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

打印：

True 
True 
orig    0.166903 
Divakar no numexpr 0.179617 
PP     0.060817 

順便說一句。如果增量t除以T可以節省更多，甚至運行完整的fft並丟棄太多。

Answer 3

這不是一個真正的另一個答案，而是@Paul Panzer的一個評論，寫作答案，因爲我需要發佈一些代碼。如果有方法可以在評論中張貼格式化好的代碼，請諮詢。

通過@保羅裝甲cumprod理念的啓發，我想出了以下內容：

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] 

雖然它似乎正確地矢量化，併產生正確的結果，它的性能更是苦不堪言。它不僅比cumprod慢得多，預計len(a)-1的指數會更多，但比原來的未引導版本要慢50％。這種糟糕表現的原因是什麼？