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