Python/numpy 切片问题解析

Question

我在使用numpy的时候遇到了一些问题。我需要一个numpy数组以一种不寻常的方式工作，也就是说，当我对数组进行切片时，它应该返回数据的视图，而不是一个副本。下面是我想做的一个例子：

假设我们有一个简单的数组，如下所示：

a = array([1, 0, 0, 0])

我想用数组中的前一个元素来更新数组中连续的条目（从左到右），使用类似这样的语法：

a[1:] = a[0:3]

这样会得到以下结果：

a = array([1, 1, 1, 1])

或者像这样：

a[1:] = 2*a[:3]
# a = [1,2,4,8]

为了进一步说明，我想要这样的行为：

for i in range(len(a)):
    if i == 0 or i+1 == len(a): continue
    a[i+1] = a[i]

不过我希望能有numpy的速度。

numpy的默认行为是对切片进行复制，所以我实际上得到的是这个：

a = array([1, 1, 0, 0])

我已经把这个数组作为ndarray的子类，所以如果需要的话，我可以对它进行进一步的修改。我只需要右侧的切片在更新左侧切片时能够持续更新。

我是在做梦吗，还是说这种魔法真的可能实现？

更新：这一切都是因为我试图使用高斯-赛德尔迭代来解决一个线性代数问题，或多或少。这是一个涉及谐函数的特殊情况，我本来不想提这个，因为这并不必要，而且可能会让事情更复杂，但我还是说出来了。

算法是这样的：

while not converged:
    for i in range(len(u[:,0])):
        for j in range(len(u[0,:])):
            # skip over boundary entries, i,j == 0 or len(u)
            u[i,j] = 0.25*(u[i-1,j] + u[i+1,j] + u[i, j-1] + u[i,j+1])

对吧？但是你可以用两种方式来做，雅可比方法是更新每个元素时不考虑已经做过的更新，直到循环结束。如果用循环的话，你会复制数组，然后从复制的数组中更新一个数组。然而，高斯-赛德尔方法使用你已经更新的信息来处理每个i-1和j-1的条目，因此不需要复制，循环应该在每次单个元素更新后“知道”这一点。也就是说，每次我们调用像u[i-1,j]或u[i,j-1]这样的条目时，之前循环中计算的信息都会存在。

我想用numpy切片来替代这种慢而复杂的嵌套循环，用一行干净的代码来实现：

u[1:-1,1:-1] = 0.25(u[:-2,1:-1] + u[2:,1:-1] + u[1:-1,:-2] + u[1:-1,2:])

但结果是雅可比迭代，因为当你切片时：u[:,-2,1:-1]你复制了数据，因此这个切片并不知道任何已做的更新。现在numpy仍然会循环，对吧？它不是并行的，只是看起来像是在python中以更快的方式循环。我想利用这种行为，试图通过某种方式“黑进”numpy，让它在我切片时返回一个指针，而不是一个副本。对吧？这样每次numpy循环时，这个切片就会“更新”，或者说只是复制更新中发生的事情。为此，我需要数组两侧的切片都是指针。

无论如何，如果有哪个聪明绝顶的人能解决这个问题，那就太棒了，但我基本上已经接受了唯一的答案就是用C语言来循环。

Answer 1

你只需要用一个循环就行。我一时间想不到有什么办法能让切片操作符像你想的那样工作，除了可能通过子类化numpy的array并用一些Python的技巧来重写相应的方法……不过更重要的是，a[1:] = a[0:3] 这个写法让我觉得完全没有道理，因为它的意思是把a的第一个值复制到接下来的三个位置。我想这会让其他人看你代码的时候感到困惑（至少前几次会这样）。

Answer 2

accumulate 是一个用来实现你想要的功能的工具；也就是说，它可以在一个数组上进行操作并把结果传播下去。下面是一个例子：

from numpy import *

a = array([1,0,0,0])
a[1:] = add.accumulate(a[0:3])
# a = [1, 1, 1, 1]

b = array([1,1,1,1])
b[1:] = multiply.accumulate(2*b[0:3])
# b = [1 2 4 8]

还有一种方法是明确指定结果数组为输入数组。下面是一个例子：

c = array([2,0,0,0])
multiply(c[:3], c[:3], c[1:])
# c = [  2   4  16 256]

Answer 3

虽然我回答得有点晚，但这个问题在谷歌上出现了，所以我想指向提问者想要的文档。你的问题很明显：在使用NumPy切片时，会创建一些临时变量。你可以用一个简单的weave.blitz调用把你的代码包裹起来，这样就能消除这些临时变量，达到你想要的效果。

想了解更多细节，可以查看PerformancePython教程中的weave.blitz部分。