numpy：評估矩陣中的函數，使用之前的數組作爲參數計算下一個

Question

我有一個m x n數組：a，其中整數m > 1E6和n <= 5。

我有功能˚F和ģ，其由這樣的：˚F（ü，ģ（ü，t））的。 ü是一個1 x n數組，t是一個標量，並且F和G返回1 x n數組。

我需要˚F評估的a每個row，並使用先前估計的行爲ü -array爲接下來的評測。我需要做出這樣的評估m。

這必須非常快。我之前對整個數組的評估印象深刻，但是這個問題需要使用先前計算的數組作爲計算下一個數組的參數。我不知道如果StringFunction可以做到這一點。

例如：

a = zeros((1000000, 4)) 
a[0] = asarray([1.,69.,3.,4.1]) 

# A is a float defined elsewhere, h is a function which accepts a float as its argument and returns an arbitrary float. h is defined elsewhere. 

def G(u, t): 
    return asarray([u[0], u[1]*A, cos(u[2]), t*h(u[3])]) 

def F(u, t): 
    return u + G(u, t) 


dt = 1E-6 

for i in range(1, 1000000): 
    a[i] = F(a[i-1], i*dt) 
    i += 1 

與上述代碼的問題是，它是緩慢的地獄。我需要用numpy毫秒來完成這些計算。

我該怎麼做我想要的？

謝謝你的時間。

親切的問候，

馬呂斯

Answer 1

這種事情在numpy中很難做到。如果我們逐列分析，我們會看到一些更簡單的解決方案。

a[:,0]是很容易的：

col0 = np.ones((1000))*2 
col0[0] = 1     #Or whatever start value. 
np.cumprod(col0, out=col0) 

np.allclose(col0, a[:1000,0]) 
True 

正如前面提到的，這將很快溢出。 a[:,1]可以沿着同樣的路線完成。

我不認爲有一種方法可以快速完成numpy內部的下兩列。我們可以求助於numba此：

from numba import auotojit 

def python_loop(start, count): 
    out = np.zeros((count), dtype=np.double) 
    out[0] = start 
    for x in xrange(count-1): 
     out[x+1] = out[x] + np.cos(out[x+1]) 
    return out 

numba_loop = autojit(python_loop) 

np.allclose(numba_loop(3,1000),a[:1000,2]) 
True 

%timeit python_loop(3,1000000) 
1 loops, best of 3: 4.14 s per loop 

%timeit numba_loop(3,1000000) 
1 loops, best of 3: 42.5 ms per loop 

雖然它的值得指出的是這個非常非常迅速收斂到pi/2並沒有多少點計算此遞歸過去〜20個值的任何初始值。這將返回完全相同的答案，雙點precision-我沒有打擾找到截止，但它是非常小於50：

%timeit tmp = np.empty((1000000)); 
     tmp[:50] = numba_loop(3,50); 
     tmp[50:] = np.pi/2 
100 loops, best of 3: 2.25 ms per loop 

你可以做的第四列類似的東西。當然，你可以autojit所有的功能，但是這給你幾個不同的選項來嘗試根據numba用法：

1. 的前兩列用cumprod
2. 對列3近似（以及可能的4），其中只有前幾個迭代計算
3. 使用autojit
4. 裹一切的autojit環（最好的選擇）
5. 的方式裏面，你已經提出了這一切行過去實施列numba 3和4〜 200將b e np.inf或np.pi/2。利用這個。

Answer 2

稍快。你的第一列基本上是2^n。計算2^n爲n高達1000000會溢出。第二列更糟糕。

def calc(arr, t0=1E-6): 
    u = arr[0] 
    dt = 1E-6 
    h = lambda x: np.random.random(1)*50.0 

    def firstColGen(uStart): 
     u = uStart 
     while True: 
      u += u 
      yield u 

    def secondColGen(uStart, A): 
     u = uStart 
     while True: 
      u += u*A 
      yield u 

    def thirdColGen(uStart): 
     u = uStart 
     while True: 
      u += np.cos(u) 
      yield u 

    def fourthColGen(uStart, h, t0, dt): 
     u = uStart 
     t = t0 
     while True: 
      u += h(u) * dt 
      t += dt 
      yield u 

    first = firstColGen(u[0]) 
    second = secondColGen(u[1], A) 
    third = thirdColGen(u[2]) 
    fourth = fourthColGen(u[3], h, t0, dt) 

    for i in xrange(1, len(arr)): 
     arr[i] = [first.next(), second.next(), third.next(), fourth.next()]