替换 NumPy 中的迭代答案

【问题标题】：Replacing iteration in NumPy替换 NumPy 中的迭代
【发布时间】：2016-02-06 16:13:42
【问题描述】：

我正在尝试重写以下例程（a、b、c、d 和 e 都是数组）：

def generate_a(b, c, e):
    a = np.zeros_like(b)
    d = np.zeros_like(b)
    for i in range(a.size):
        a[i] = (b[i] / c[i]) + d[i-1]
        d[i] = a[i] * e[i]
    return a

不使用“for”循环，因为例程需要运行数百万次。这是一个小技巧，因为 d 中任何给定单元格的值取决于为 a 中的单元格计算的结果，而结果又取决于为 d 计算的最后一个值。关于如何在没有迭代的情况下完成此操作的任何想法？

【问题讨论】：

对于元素a[0]，预计使用d的什么值？

标签： python arrays numpy recursion vectorization

【解决方案1】：

受solution 的启发，这篇文章表明递归是vectorizable，如果递归可以跟踪。正如链接解决方案中所述，这不应该很快到目前为止，但可以利用并行处理器，这种方法可以尝试。因此，请将这篇文章作为有关如何跟踪和矢量化递归的指导。这是实现 -

def generate_a_vectorized(b, c, e):

    K = (b/c)*e
    N = e.size

    mask = np.linspace(N,1,N,dtype=int)[:,None] <= np.arange(N)

    Pn = np.tile(e[None],(N,1))
    Pn[mask] = 1
    En = Pn[:,::-1].cumprod(1)[:,::-1]
    En[mask] = 0
    An = np.append(0,K[:-1])

    d_out = En.dot(An)[::-1] + K
    return (b/c) + np.append(0,d_out[:-1])

示例运行和验证输出 -

In [279]: M = 50
     ...: b = np.random.rand(M)
     ...: c = np.random.rand(M)
     ...: e = np.random.rand(M)
     ...: 

In [280]: np.allclose(generate_a(b, c, e),generate_a_vectorized(b, c, e))
Out[280]: True

In [281]: %timeit generate_a(b,c,e)
10000 loops, best of 3: 79.4 µs per loop

In [282]: %timeit generate_a_vectorized(b,c,e)
10000 loops, best of 3: 157 µs per loop

【讨论】：

【解决方案2】：

因为a[i] 依赖于d[i-1] 而d[i-1] 依赖于a[i-1]，所以a[i] 也依赖于a[0], a[1], ... a[i-1]，因此您将无法按原样并行化计算。您需要接受该循环或找出一种封闭形式的方法来进行计算。

【讨论】：

封闭形式是可计算的

【解决方案3】：

您可以使用 map 和 lambda。

试试：

import numpy as np

def func(i,l):
    l[0][i] = (l[1][i] / l[2][i]) + l[3][i-1]
    l[3][i] = l[0][i] * l[4][i]

    return l[0][i]

def generate_a(b, c, e):
    a = np.zeros_like(b)
    d = np.zeros_like(b)

    a = map(lambda i:func(i,[a, b, c, d, e]), xrange(a.size))

    return a

print generate_a([5, 6, 7],[1,2,3],[8,9,12])

【讨论】：

这样更优雅，但是由于 map 函数仍然是迭代的，所以并没有节省多少时间。

【解决方案4】：

一个关于 Mathematica 的快速想法（因为b[i]/c[i] 独立于计算，我在下面将其称为bc[i]）：

Wolfram Language (Raspberry Pi Pilot Release)
Copyright 1988-2015 Wolfram Research
Information & help: wolfram.com/raspi

In[1]:= RSolve[ { a[i]==bc[i] + a[i-1]*e[i-1], a[-1]==0 }, a[i], i]                    

Out[1]= {{a[i] -> 

>      Product[e[K[1]], {K[1], 1, -1 + i}] 

                      bc[1 + K[2]]
>       Sum[---------------------------------, {K[2], -1, -1 + i}]}}
            Product[e[K[1]], {K[1], 1, K[2]}]

正如 Randy C 之前所说的那样，除非之前的计算会产生更简单的结果，否则我强烈怀疑你能否得到比你已经拥有的更好的结果。 你可以获得更好的代码，但可能没有更高效的代码。

【讨论】：