python - numpy矩阵代数最佳实践

Question

我的问题是关于下面最后一行：mu@sigma@mu。为什么有效？一维nDarray被视为行向量还是列向量？不管怎样，不是应该是mu.T@sigma@mu还是mu@sigma@mu.T？我知道mu.T仍然返回mu，因为mu只有一个维度，但是解释器似乎还是太聪明了。

>> import numpy as np
>> mu = np.array([1, 1])
>> print(mu)

[1 1]

>> sigma = np.eye(2) * 3
>> print(sigma)

[[ 3.  0.]
 [ 0.  3.]]

>> mu@sigma@mu

6.0

更一般地说，这是python中矩阵代数的一个更好的实践：使用ndarray和 @来执行上面的矩阵乘法（更干净的代码），或者使用 np.matrix和下面的重载 *（数学上不那么混乱）

>> import numpy as np
>> mu = np.matrix(np.array([1, 1]))
>> print(mu)

[[1 1]]

>> sigma = np.matrix(np.eye(2) * 3)
>> print(sigma)

[[ 3.  0.]
 [ 0.  3.]]

>> a = mu * sigma * mu.T
>> a.item((0, 0))

6.0

Answer 1

python执行从左到右的链接操作：

In [32]: mu=np.array([1,1])
In [33]: sigma= np.array([[3,0],[0,3]])
In [34]: mu@sigma@mu
Out[34]: 6

与执行两个表达式相同：

In [35]: temp=mu@sigma
In [36]: temp.shape
Out[36]: (2,)
In [37]: temp@mu
Out[37]: 6

在我的评论（删除）中，我声称 @只是在做 np.dot。那不太对。文档以不同的方式描述1d数组的处理。但最终的形状是一样的：

In [38]: mu.dot(sigma).dot(mu)
Out[38]: 6
In [39]: mu.dot(sigma).shape
Out[39]: (2,)

对于一维和二维阵列， np.dot和 @应该产生相同的结果它们在处理高维数组方面有所不同。
从历史上看，numpy使用数组，可以是0d、1d和on-up。 np.dot是最初的矩阵乘法方法/函数。
np.matrix被添加，主要是为了方便任性的MATLAB程序员它只允许二维数组（就像旧的，90年代的MATLAB一样）它超载 __mat__（*）

def __mul__(self, other):
    if isinstance(other, (N.ndarray, list, tuple)) :
        # This promotes 1-D vectors to row vectors
        return N.dot(self, asmatrix(other))
    if isscalar(other) or not hasattr(other, '__rmul__') :
        return N.dot(self, other)
    return NotImplemented

Mu*sigma和 Mu@sigma的行为相同，尽管调用树不同

In [48]: Mu@sigma@Mu
...
ValueError: shapes (1,2) and (1,2) not aligned: 2 (dim 1) != 1 (dim 0)

Mu*sigma产生（1,2）矩阵，该矩阵不能与（1,2）相乘，因此需要转置：

In [49]: Mu@sigma@Mu.T
Out[49]: matrix([[6]])

注意这是一个（1，1）矩阵。如果需要标量，则必须使用 item。（在Matlab中没有标量这样的东西。所有东西都有形状/大小。）
@是python和numpy的一个相对较新的添加。它是作为一个未实现的操作符添加到python中的。numpy（可能还有其他包）已经实现了它。
它使链式表达式成为可能，尽管我对 dot中的链式 [38]没有任何问题。它在处理高维情况时更有用。
这意味着使用旧的 np.matrix类的理由更少了。（类似矩阵的行为在 scipy.sparse矩阵类中根深蒂固。）
如果你想要“数学纯度”，我建议你采用数学物理的方法，并使用爱因斯坦符号-如 np.einsum中所实现的。
在数组这么小的情况下，计时反映的调用结构比实际的计算数量还要多：

In [57]: timeit mu.dot(sigma).dot(mu)
2.79 µs ± 7.75 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
In [58]: timeit mu@sigma@mu
6.29 µs ± 31.4 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
In [59]: timeit Mu@sigma@Mu.T
17.1 µs ± 134 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
In [60]: timeit Mu*sigma*Mu.T
17.7 µs ± 517 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

请注意，“老式的” dot是最快的，而两个矩阵版本都比较慢。
论转置
我可以补充说，在一个类似于 mu@sigma@mu.T的表达式中，首先计算 .T，因此 @（ matmul）不知道有人试图转置 mu。它只是使用转置的结果。记住，解析表达式的是python解释器，而不是 numpy函数。因此，如果 mu.T对 mu没有任何作用，就像一维数组一样， .T也不会对 @产生影响。
transpose对于二维数组定义良好。 numpy转置是以这样的方式编写的，它也可以用于任何维度的数组它改变轴的顺序，但从不添加维度还有其他的工具，比如 reshape和 np.newaxis索引。
在八度音阶中，转置仅为二维对象定义

>> ones(2,3,4)'
error: transpose not defined for N-d objects

一维对象不存在于Matlab中。

In [270]: np.ones((2,3,4),int).T.shape
Out[270]: (4, 3, 2)

内外产品
np.dot非常明确地说明了它如何处理1d数组，“向量”：
如果 a和 b都是一维数组，则它是向量的内积。
（没有复杂的共轭）。
在 matmul中，描述更加复杂，但结果相同。
关于你的评论：
a=[1，2]和b=[3，5]为（2，）ndarray，a@b可以表示[1，2]*[3，5]'=13，也可以表示[1，2]*[3，5]=[[3，5]，[6，10]]
在 numpy中使用产品的各种方法有：

In [273]: A = np.array([1,2]); B = np.array([3,5])

元素乘法（在Matlab中 .*）

In [274]: A*B
Out[274]: array([ 3, 10])

内积

In [275]: A@B       # same as np.dot(A,B)
Out[275]: 13

外部产品

In [276]: np.outer(A,B)
Out[276]: 
array([[ 3,  5],
       [ 6, 10]])

另一种获得内部产品的方法是：

In [278]: np.sum(A*B)
Out[278]: 13

用爱因斯坦符号（来自数学物理）表示内部和外部：

In [280]: np.einsum('i,i',A,B)
Out[280]: array(13)
In [281]: np.einsum('i,j',A,B)
Out[281]: 
array([[ 3,  5],
       [ 6, 10]])

利用广播获取外部

In [282]: A[:,None]*B[None,:]
Out[282]: 
array([[ 3,  5],
       [ 6, 10]])

您所期望的 [1, 2]'在 numpy中实现为 A[:,None]，将1d数组转换为列向量（2d）。
八度音阶在 numpy之后很长一段时间增加了广播我不知道MATLAB有没有那么先进：）
@不进行广播，但可以使用列向量或行向量：

In [283]: A@B[:,None]       # your imagined A*B'
Out[283]: array([13])

要使用 @获取外部，我需要添加一个维度，使 A成为列向量，并使 B成为行向量：

In [284]: A[:,None]@B
ValueError: shapes (2,1) and (2,) not aligned: 1 (dim 1) != 2 (dim 0)
In [285]: A[:,None]@B[None,:]
Out[285]: 
array([[ 3,  5],
       [ 6, 10]])

当数组真的是二维的，包括一维为1的情况时， @/ dot的行为与常见的矩阵约定没有太大区别当维数为1或大于2时，来自MATLAB的直觉就失败了，部分原因是MATLAB不能处理这些问题。
这个倍频程误差表明，n-d矩阵只依赖于2d矩阵。它们不是真正的n-d。

>> ones(2,3,4) * ones(2,3,4)
error: operator *: nonconformant arguments (op1 is 2x12, op2 is 2x12)