Numpy 计算协方差的结果应该如何理解？

[Richard14]

已知两向量 a 和 b ，利用协方差公式计算其相关性

import numpy as np 
a = np.array([3,5,4,12,9])
b = np.array([5,15,5,6,7])

根据公式查得协方差公式为 Cov(X,Y)=E(XY)-E(X)E(Y)

那么E(X)=(3+5+4+12+9)/5=6.6

E(Y)=(5+15+5+6+7)/5=7.6

E(XY)=(3*5+5*15+4*5+12*6+9*7)/5=49

协方差结果应该等于 -1.16

但是 numpy 中给出的结算结果是

>>> np.array([a,b])
[[14.3  -1.45]
 [-1.45 17.8 ]]

为什么协方差计算结果不是一个值而是 2x2 的矩阵呢？而且其中也没有结果等于-1.16 ，应该如何理解 numpy 的计算结果？

另外需求上是提供一个图片指纹（一维向量），想要尝试使用协方差计算最近似图片的思路，也就是和数据库里备选的几万个同样长度的向量分别做协方差，并筛选出相似度最高的。这种情况下应该如何写 numpy 会比较快呢？

[deanguqiang]

1. 你用 numpy 计算的是协方差矩阵
2. 协方差矩阵 C12/C21 的就是你说的协方差
3. 为什么协方差和你手算的不一样？我估计是 Numpy 用 N-1 来平均的，-1.16/4*5=-1.45
4. 你真正需要的可能是相关系数

[noqwerty]

np.cov([a, b], bias=True)

[jaredyam]

关于#1 的 4：

协方差矩阵也可以作为图片的特征矩阵，从统计角度计算图片相似性。前提是你需要假设它的概率分布，从而根据最大似然估计推导依赖协方差矩阵的相似性度量。

[Richard14]

@deanguqiang @noqwerty 是的，相关系数也准备尝试，但是因为相关系数基于协方差，而协方差计算结果不准确，所以先发帖问了一下。将 bias 设置为 true 后计算结果正常，但是根据我的理解输出矩阵的 m 行 n 列的意思应该代表第 m 行的向量与第 n 行的向量的协方差，计算复杂度是向量总个数的平方 /2 ，当输入量大时（比如十万个向量）感觉浪费了很多算力（实际只需要衡量唯一输入与所有备选的相似度，输出结果只需要取当前结果的第一行就可以了），这种计算应该怎么写呢？

@jaredyam 是的，在网上查找资料时发现协方差似乎有两种公式，上文述是直接用期望计算，还有一种需要考虑到数据的分布概率，但是网上相关资料似乎大多是不考虑概率的计算公式。想要加入概率的话应该怎么设计计算呢，能够 numpy 化吗？

[necomancer]

np.cov(x, y, ddof=0)
# 11.44 1.16
# 1.16 14.24

cov 是所有的都算的，矩阵是 xx, xy, yx, yy ，你需要的仅是 xy ，你就直接 np.mean((a-a.mean())*(b-b.mean())) 就行了。另外，关联性一般用 correlate ，协方差(ddof=0)应该是关联函数的第一个值，后面可以看到衰减。

[Richard14]

@necomancer 感谢回复，如果我需要判断 10 个向量里与输入向量最接近的，除了使用 for 循环挨个计算以外，有什么快捷写法吗

[necomancer]

@Richard14 假设十个向量的形状是 a=(10, N)，目标向量是 b=(N,)，你可以用 np.einsum：

ret_index = np.argmin(np.einsum('ij,j->i', a,b)/np.linalg.norm(a, axis=-1)/np.linalg.norm(b, axis=-1))
ret_vec = a[ret_index]

这个简洁但是会比较慢，可以考虑用 numba 的 guvectorize
from numba import float64, guvectorize
@guvectorize([(float64[:, :], float64[:], float64[:])],'(n,p),(p)->(n)', target='parallel') # target = 'cpu', 'gpu'
....def batch_dot(a, b, ret):
........for i in range(a.shape[0]):
........t1=t2=t3=0
............for j in range(a.shape[1]):
................t1 += a[i, j] * b[j]
................t2 += a[i, j] * a[i, j]
................t3 += b[j] * b[j]
............ret[i] = t1 / t2**0.5/t3**0.5
ret = batch_dot(a,b)

[Richard14]

@necomancer 感叹于 numpy 的深度

[necomancer]

@Richard14 我的回复好像有点脑残……你是要判断相似性或者说距离的话算法是不一样的，如果是距离那么直接 np.argmin(np.sum((a-b) ** p, axis=-1))，比如 p=2 是欧式距离。我写的是 cos(theta)，做距离的话一般是 1 - cos(theta) 啥的，这个你自己调整一下吧。numpy 是很炫酷的。