在Python中，矩阵的就地赋值是什么意思？

Question

假设我按照以下方式初始化一个矩阵：

import scipy
m = scipy.zeros((10, 10))

现在我需要做一些计算，并将结果赋值给m。在这个赋值过程中，m的大小不会改变，因此我认为如果能够直接在原地完成赋值，速度会更快。

m = scipy.array([[i * j for j in range(10)] for i in range(10)])

我担心在上面的代码中，会创建一个临时矩阵来保存结果，然后将m赋值为这个值。这是低效的，因为它涉及到分配一个新的矩阵。更有效的解决方案是直接将值存储在m中，可以这样表达：

for i in range(10):
    for j in range(10):
        m[i,j] = i * j

但是假设生成器表达式对我的代码编写方式更加方便，我想知道的是：在上述生成器表达式中，我是否进行了额外的矩阵分配？

- a06e

---

数组分配可能看起来很昂贵，如果您习惯于像C这样的语言，但是与CPython的无JIT字节码解释和动态执行相比，它的成本非常低廉。如果您希望NumPy/SciPy代码高效运行，最关键的是使用NumPy向量化操作将工作推入避免所有开销的C级别循环中。首先尝试最小化分配不会有太大帮助，甚至可能使您的代码变慢。 - user2357112

Answer 1

让我们进行一些实际的时间测试：

In [793]: timeit m=np.array([[i*j for j in range(N)] for i in range(M)])
10000 loops, best of 3: 47.8 µs per loop
In [794]: %%timeit
   .....: m=np.zeros((N,M),int)
   .....: for i in range(M):
    for j in range(N):
        m[i,j] = i*j
   .....: 
10000 loops, best of 3: 40.2 µs per loop

因此，预分配和赋值速度略快 - 但并不是非常明显。

与向量乘法相比:

In [796]: timeit np.arange(M)[:,None]*np.arange(N)[None,:]
10000 loops, best of 3: 17.1 µs per loop

对于更大的数组，执行相同的操作：

In [797]: N,M=1000,1000
In [798]: timeit m=np.array([[i*j for j in range(N)] for i in range(M)])
1 loops, best of 3: 325 ms per loop
In [799]: %%timeit
m=np.zeros((N,M),int)
for i in range(M):
    for j in range(N):
        m[i,j] = i*j
   .....: 
1 loops, best of 3: 338 ms per loop
In [800]: timeit np.arange(M)[:,None]*np.arange(N)[None,:]
100 loops, best of 3: 12.5 ms per loop

这两次迭代保持了领先优势；向量化更好。

我可以通过使用 fromiter 来节省一些迭代的时间，但是远不及向量化。

In [805]: timeit np.fromiter([i*j for j in range(N) for i in range(M)],int).reshape(N,M)
1 loops, best of 3: 235 ms per loop

这是一个常见的问题，我太懒了，不想去寻找最佳的重复解决方案。:) 通常人们声称他们的计算是一些只接受标量的复杂黑匣子，所以没有办法对其进行矢量化。

有一个np.vectorize函数可以包装你的计算，但它旨在简化诸如广播之类的事情，并不保证加速代码。它仍然需要迭代。

如果计算量很小且速度较快，值得注意迭代方法；但如果计算很复杂，则花费在迭代机制上的时间比例很小，应该关注黑匣子的速度。

- hpaulj

Answer 2

你第一个解决方案，列表推导式的问题在于它生成了一个嵌套列表并将其分配给m。然而，从你的第一条语句开始，你似乎希望m成为一个numpy数组（这是在执行scipy.zeros()时创建的内容）。因此，你实际上创建了一个数组，然后用一个列表覆盖了它。如果你想保持数据结构为np.array，则你的嵌套for循环是最好的方法。

另外，你说“matrix”，但创建了一个数组。如果你想要一个真正的矩阵（例如进行矩阵运算），请将你的嵌套列表推导式传递给np.matrix()：

# assuming you've already run `import numpy as np`
In [5]: m = np.matrix([[i * j for j in range(10)] for i in range(10)])

In [6]: m
Out[6]: 
matrix([[ 0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
        [ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9],
        [ 0,  2,  4,  6,  8, 10, 12, 14, 16, 18],
        [ 0,  3,  6,  9, 12, 15, 18, 21, 24, 27],
        [ 0,  4,  8, 12, 16, 20, 24, 28, 32, 36],
        [ 0,  5, 10, 15, 20, 25, 30, 35, 40, 45],
        [ 0,  6, 12, 18, 24, 30, 36, 42, 48, 54],
        [ 0,  7, 14, 21, 28, 35, 42, 49, 56, 63],
        [ 0,  8, 16, 24, 32, 40, 48, 56, 64, 72],
        [ 0,  9, 18, 27, 36, 45, 54, 63, 72, 81]])

即使你最终需要一个数组，也可以像上面一样将嵌套的列表解析器传递给数组构造函数，然后就可以了。

- MattDMo

---

抱歉，我在赋值中漏掉了scipy.array(...)。我已经编辑了问题。那么你是说赋值比for循环更好？哪个更快？ - a06e

我不建议使用np.matrix - 除非OP执着于旧的MATLAB操作。 - hpaulj

@hpaulj 为什么这样？ - MattDMo

np.matrix 是 np.ndarray 的子类，它始终是二维的。这使得某些事情更简单，但也有一定的限制。它的矩阵乘法表达式更简单，但逐元素乘法更复杂。 - hpaulj

最近的一个例子是http://stackoverflow.com/questions/35117660/numpy-ndarray-multiplication-switching-to-matrix-multiplication，在这里无意中切换到“matrix”会导致问题。 - hpaulj

@hpaulj 很有趣，感谢您的解释。所以基本上带走的教训是尽可能不要使用 np.matrix 吗？ - MattDMo

Answer 3

第二个任务（生成器）确实创建了一个新的矩阵。 如果您使用Python的id()函数，您会发现在这个任务之后m指向了一个不同的位置。

例如：

>> import scipy
>> m = scipy.zeros((10, 10))
>> id(m)
4455211696
>> m = scipy.array([[i * j for j in range(10)] for i in range(10)])
>> id(m)
4478936688

- Tomer Levinboim