我的问题是关于下面最后一行:mu@sigma@mu
。为什么有效?一维nDarray被视为行向量还是列向量?不管怎样,不是应该是mu.T@sigma@mu
还是mu@[email protected]
?我知道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
最佳答案
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@[email protected]
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@[email protected]
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@[email protected]
的表达式中,首先计算.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)
关于python - numpy矩阵代数最佳实践,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/48492429/