numpy中的最优（广播）矩阵除法。是否避免临时数组？

Question

Numpy允许不同大小的矩阵进行加法、乘法和除法，只要遵循一些特定的广播规则。另外，创建临时数组会大大影响Numpy的速度。

下面的timing结果让我感到惊讶……这是怎么回事呢？

In [41]: def f_no_dot(mat,arr):
   ....:     return mat/arr

In [42]: def f_dot(mat,arr):
   ....:     denominator = scipy.dot(arr, scipy.ones((1,2)))
   ....:     return mat/denominator

In [43]: mat = scipy.rand(360000,2)

In [44]: arr = scipy.rand(360000,1)

In [45]: timeit temp = f_no_dot(mat,arr)
10 loops, best of 3: 44.7 ms per loop

In [46]: timeit temp = f_dot(mat,arr)
100 loops, best of 3: 10.1 ms per loop

我原以为f_dot会更慢，因为它需要创建一个临时数组作为分母，而我认为f_no_dot这个步骤是省略掉的。我还要提到，对于f_no_dot，这些时间是线性增长的（随着数组大小，最多到10亿），而f_dot的增长速度稍微差于线性。

Answer 1

我原以为 f_dot 会比较慢，因为它需要创建一个临时数组 denominator，而我假设 f_no_dot 这个步骤是跳过的。

其实，创建临时数组确实是被跳过了，这就是为什么 f_no_dot 会比较慢（但占用的内存更少）。

对同样大小的数组进行逐元素操作会更快，因为 numpy 不用担心数组的步幅（维度、大小等）。

使用广播的操作通常会比不使用广播的操作稍微慢一些。

如果你的内存足够，创建一个临时副本可以让你速度更快，但会占用更多内存。

比如，比较这三个函数：

import numpy as np
import timeit

def f_no_dot(x, y):
    return x / y

def f_dot(x, y):
    denom = np.dot(y, np.ones((1,2)))
    return x / denom

def f_in_place(x, y):
    x /= y
    return x

num = 3600000
x = np.ones((num, 2))
y = np.ones((num, 1))


for func in ['f_dot', 'f_no_dot', 'f_in_place']:
    t = timeit.timeit('%s(x,y)' % func, number=100,
            setup='from __main__ import x,y,f_dot, f_no_dot, f_in_place')
    print func, 'time...'
    print t / 100.0

这和你的结果时间差不多：

f_dot time...
0.184361531734
f_no_dot time...
0.619203259945
f_in_place time...
0.585789341927

不过，如果我们比较内存使用情况，就会更清楚了……

我们的 x 和 y 数组的总大小大约是 27.5 + 55 MB，也就是 82 MB（对于 64 位整数）。另外还有大约 11 MB 的额外开销，比如导入 numpy 等。

返回 x / y 作为一个新数组（也就是说，不是做 x /= y）会需要另一个 55 MB 的数组。

运行 f_dot 100 次： 我们在这里创建了一个临时数组，所以我们预计内存使用量会是 11 + 82 + 55 + 55 MB，约 203 MB。结果确实是这样……

运行 f_no_dot 100 次： 如果没有创建临时数组，我们预计峰值内存使用量会是 11 + 82 + 55 MB，约 148 MB……
……这正是我们看到的结果。

所以，x / y 并没有创建额外的 num x 2 临时数组来进行除法运算。

因此，除法的速度比在两个同样大小的数组上操作要慢很多。

运行 f_in_place 100 次： 如果我们可以就地修改 x，那么如果内存是主要考虑因素，我们可以节省更多内存。

基本上，numpy 在某些情况下会为了节省内存而牺牲速度。

Answer 2

你看到的情况很可能是因为在处理小的 (2,) 维度时，迭代的开销比较大。Numpy（1.6版本之前）在处理只有相邻数组（形状相同的数组）时效率比较高。随着最后一个维度的大小增加，这种影响会减小。

要了解相邻性（contiguity）的影响，可以看这个例子：

In [1]: import numpy
In [2]: numpy.__version__
Out[2]: '1.5.1'
In [3]: arr_cont1 = numpy.random.rand(360000, 2)
In [4]: arr_cont2 = numpy.random.rand(360000, 2)
In [5]: arr_noncont = numpy.random.rand(360000, 4)[:,::2]
In [6]: arr_bcast = numpy.random.rand(360000, 1)
In [7]: %timeit arr_cont1 / arr_cont2
100 loops, best of 3: 5.75 ms per loop
In [8]: %timeit arr_noncont / arr_cont2
10 loops, best of 3: 54.4 ms per loop
In [9]: %timeit arr_bcast / arr_cont2
10 loops, best of 3: 55.2 ms per loop

。不过在Numpy 1.6.0及以上版本中，这种情况有了很大的改善：

In [1]: import numpy
In [2]: numpy.__version__
Out[2]: '1.6.1'
In [3]: arr_cont1 = numpy.random.rand(360000, 2)
In [4]: arr_cont2 = numpy.random.rand(360000, 2)
In [5]: arr_noncont = numpy.random.rand(360000, 4)[:,::2]
In [6]: arr_bcast = numpy.random.rand(360000, 1)
In [7]: %timeit arr_cont1 / arr_cont2
100 loops, best of 3: 5.37 ms per loop
In [8]: %timeit arr_noncont / arr_cont2
100 loops, best of 3: 6.12 ms per loop
In [9]: %timeit arr_bcast / arr_cont2
100 loops, best of 3: 7.81 ms per loop

。（上面所有的时间测量可能只准确到1毫秒。）

另外要注意，临时变量的开销并没有那么大：

In [82]: %timeit arr_cont1.copy()
1000 loops, best of 3: 778 us per loop

编辑：上面提到的 arr_noncont 其实也是有点相邻的，步长是 2*itemsize，所以内部循环可以被展开——Numpy可以让它的速度接近相邻数组的速度。如果使用广播（broadcasting）或者像 numpy.random.rand(360000*2, 2)[::2,:] 这样非常不相邻的数组，内部循环就无法展开，这些情况会稍微慢一些。如果Numpy能为每个循环即时生成专门的机器代码，那改善这个问题是有可能的，但目前它还没有做到（至少现在还没有 :)