混合使用numpy矩阵和数组的危险性

Question

我正在处理一个与科学/工程相关的应用程序，其中有很多线性代数矩阵乘法，因此我使用Numpy矩阵。然而，在Python中有许多函数可以交替接受矩阵或数组类型。听起来不错，是吧？实际上并不是这样。让我通过一个例子来说明问题：

from scipy.linalg import expm
from numpy import matrix

# Setup input variable as matrix
A = matrix([[ 0, -1.0,  0,  0],
            [ 0,  0,  0,  1.0],
            [ 0,  0,  0,  0],
            [ 0,  0,  1.0,  0]])

# Do some computation with that input
B = expm(A)

b1 = B[0:2, 2:4]
b2 = B[2:4, 2:4].T

# Compute and Print the desired output
print "The innocent but wrong answer:"
print b2 * b1

print "The answer I should get:"
print matrix(b2) * matrix(b1)

当你运行时，会得到以下结果：

The innocent but wrong answer:
[[-0.16666667 -0.5       ]
 [ 0.          1.        ]]
The answer I should get, since I expected everything to still be matrices:
[[ 0.33333333  0.5       ]
 [ 0.5         1.        ]]

如何避免这种混淆？在变量周围不断地添加matrix()调用来确保它们仍然是矩阵会让代码变得非常混乱。似乎在这方面没有标准，因此可能会导致难以检测的错误。

- Hamid

---

这就是为什么人们使用Java和其他静态语言的原因。IDE和编译器会因为使用不同类型而让你头疼不已（而IDE会告诉你类型）。 - Snakes and Coffee

如果你想要点积，最好使用numpy.dot而不是依赖于矩阵重载乘法运算符。明确比隐含更好。 - David Cain

在像这样的方程式中使用np.dot真的很烦人，例如：A = GPG_tran+ MUM_tran。 - Mehdi

@Mehdi，我同意你的看法，但是我发现坚持使用numpy数组总体上是更好的体验。不用担心返回类型，而且在矩阵上使用squeeze来减少维度（例如[[1]]）可能会很麻烦。 - Hamid

在今年的一次会议上，他们提到numpy.matrix将很快被淘汰，因为一种替代稀疏表示方法已经被采用。这是避免使用numpy.matrix的另一个原因。 - Hamid

Answer 1

我倾向于在numpy中使用array而不是matrix，原因如下：

1. matrix严格限制为2D，而您可以拥有任意维度的numpyarray。
2. 除了一些差异外，对于Matlab用户来说，array和matrix操作基本上是interchangeable的。
3. 如果您始终使用array，则会使用numpy.dot()（或Python 3.5中的新@二进制运算符）进行矩阵乘法。这将防止在代码中不确定*实际执行的问题。当遇到乘法错误时，您可以更轻松地找到问题，因为您确定正在尝试执行什么类型的乘法。

因此，我建议您尽量坚持使用numpy.array，但也要记住array和matrix之间的区别。

最后，我发现在bpython上使用numpy/scipy非常愉快。自动提示帮助你更快地学习要使用的函数的属性，而不必不断查阅numpy/scipy文档。 编辑： array与matrix之间的区别可能最好在此处回答："'array'或'matrix'？我应该使用哪个？"

- K Z

---

但是函数的可读性会受到牺牲。例如： FGQG.TF.T**2 比任何点(a,b)等价物更易读。你不同意吗？此外，如果你还没有注意到，我是一个Matlab转换者，并且认为A * B被定义为矩阵乘法非常自然。 - Hamid

@user1609675 我认为这是一个很有道理的观点。但由于我之前所有的科学计算都是用Matlab完成的，所以语法差异对我来说几乎可以忽略不计，因为我必须重新学习语法。需要一段时间才能适应“dot()”，但之后我就不再被它困扰了。它变得非常自然 :) - K Z

@user1609675 是的，我曾经在机器学习和自然语言处理方面使用过大量的Matlab。我认为重新学习每个语法并逐渐忘记Matlab的方式，并通过重构心理模型来改变思维方式是有帮助的。这需要一段时间，但我认为这是值得的。 - K Z

你可以使用np.einsum来干净地完成FGQG.TF.T**2。虽然这意味着放弃使用特定优化的矩阵乘法库，但也意味着放弃了很多临时变量。具体情况而定，可能会更快。 - Eelco Hoogendoorn

@EelcoHoogendoorn 这会是什么样子呢？np.einsum对我来说是新的。 - Hamid

我不确定我能比文档更好地解释它；它使用爱因斯坦符号评估任意张量积，示例代码在此：http://docs.scipy.org/doc/numpy/reference/generated/numpy.einsum.html三个二阶张量的乘积，等价于多次np.dot操作，看起来像这样：np.einsum('ij,jk,kl->il', a,b,c)，其中重复出现的j和k索引是“缩并”的（求和）。它允许您以高效且易读的方式表达许多矩阵运算；但是否更有效取决于具体情况。但它无疑更加灵活。 - Eelco Hoogendoorn

Answer 2

混合矩阵和常规ndarrays确实可能很棘手，而且通常不值得麻烦。我赞同其他帖子的建议，建议您坚持使用数组。

尽管如此，在您的特定示例中，问题来自于expm。根据文档，它以常规ndarray作为参数并输出一个ndarray。如果您想将输出转换回matrix，可以使用：

B = matrix(expm(A))

或者

B = expm(A).view(matrix)

现在，B是一个矩阵，B的切片本身也是矩阵，您的乘法将按预期工作。

因此，建议始终检查函数输出的类型。

- Pierre GM