如何对要应用于numpy矩阵的每一行的函数进行numpy矢量化

Question

我写了这个函数，用于计算给定多变量高斯参数的x的概率。其中，x是具有2个特征的特征，μ是具有两个特征的向量，西格玛是2x2。

def prob(x, mu, sigma):
n = len(x)
var = x - mu
sigma_inv = np.linalg.inv(sigma)
rhs = np.exp(-0.5*np.matmul(np.matmul(var.T, sigma_inv),var))
lhs = 1/(((2*np.pi)**(n/2))*np.linalg.det(sigma)**.5)
return rhs*lhs

但这只适用于x是1D数组的情况。我希望能够对我目前拥有的多维x(例如x是100x2(进行矢量化和优化。

for i in range(len(x)):
curr_prob = prob(x[i], mu, sigma)

if i == 0:
prob = curr_prob
else:
prob = np.append(prob, curr_prob)

但这是非常缓慢的。我听说有一种方法可以使用np.vectorize或np.pyfunc，但我不知道如何应用它们。

Answer 1

使用numpy.apply_along_axis:怎么样

np.apply_along_axis(prob, 1, x, mu, sigma)

Answer 2

查看形状如何在当前函数中流动

x (2,), mu (2,), sigma (2,2)
def prob(x, mu, sigma):
n = len(x)             # 2
var = x - mu           # (2,)-(2,)
sigma_inv = np.linalg.inv(sigma)    # (2,2) no change w/ x
rhs = np.exp(-0.5*np.matmul(np.matmul(var.T, sigma_inv),var))
# var is (2,), var.T is the same (2,)
# (2,)@(2,2)=>(2,); (2,)@(2,)=> scalar
# np.exp(-0.5 * var @ sigma_inv @ var
lhs = 1/(((2*np.pi)**(n/2))*np.linalg.det(sigma)**.5)
# np.linalg.def(sigma)**.5 - no dependence on x
# the whole lhs doesn't vary with x
# lhs is scalar
return rhs*lhs    # scalar

现在考虑如果x是(100,2(会发生什么变化

def prob(x, mu, sigma):
n = x.shape[-1]    
sigma_inv = np.linalg.inv(sigma)    # (2,2) or (n,n)?
lhs = 1/(((2*np.pi)**(n/2))*np.linalg.det(sigma)**.5) # scalar
var = x - mu           # (100,2)-(2,)
# by broadcasting this is (100,2)-(1,2)=>(100,2) 
# no change needed
rhs = var@sigma_inv@var
# (100,2) @ (2,2) => (100,2)
# (100,2) @ (100,2)  oops
# var@sigma_inv@var.T   
# (100,2) with (2,100)=>(100,100)   no! 
# np.einsum('ij,jk,ik->i',var,sigma_inv,var) 
# 
rhs = np.exp(-0.5*rhs)    # (100,2)
return rhs*lhs

在我的第一次尝试中，var@sigma_inv@var与(2，(一起工作，但与(n，2(一起出现错误。但是einsum表达式得到了i，批处理维度。我可能也可以纠正双@来纠正这个错误。

def prob(x, mu, sigma):
n = x.shape[-1]    
sigma_inv = np.linalg.inv(sigma)
lhs = 1/(((2*np.pi)**(n/2))*np.linalg.det(sigma)**.5)
var = x - mu
rhs = np.einsum('ij,jk,ik->i',var,sigma_inv,var)
rhs = np.exp(-0.5*rhs)
return rhs*lhs

测试：

In [495]: x = np.array([1,2.]); mu=np.array([.5,.5]);
In [496]: sigma = np.array([[1,.5],[.5,3]])
In [497]: X = np.array([x,x,x+1])
In [498]: prob(X,mu,sigma)
Out[498]: array([0.06375113, 0.06375113, 0.01785457])

您的函数还为x生成0.06375。

einsum可以写成：

np.squeeze(var[:,None,:]@sigma_inv@var[:,:,None])

但这并不清楚，而且可能在速度上也没有太大区别。