python - 向量化稀疏矩阵的 hellinger - NumPy/Python

Question

我正在寻找 Hellinger distance在单个分布 p 和稀疏矩阵 dist_mat 的每一行之间。我想返回一个维度为 1*N 的向量，其中 N 是 dist_mat 中的行数。

def hellinger(p, dist_mat):
    return np.sqrt(1/2) * np.sqrt(  np.sum((np.sqrt(p) - np.sqrt(dist_mat))**2)  )

使用上面的函数，如果我们尝试一个测试用例:

row = np.array([0, 0, 1, 2, 2, 2])
col = np.array([0, 2, 2, 0, 1, 2])
data = np.array([1, 2, 3, 4, 5, 6])
csr_matrix((data, (row, col)), shape=(3, 3)).toarray()
test = np.array([0,21,0])
hellinger(test,csr_matrix((data, (row, col)), shape=(3, 3)))
>>> 4.3633103660024926

它返回一个标量，而不是一个向量。因此，对于上面的示例，我想要一个包含 hellinger 距离的结果列表。像这样的东西:

hellinger(test,csr_matrix((data, (row, col)), shape=(3, 3)))
>>> [3.46,3.46,2.78] # hellinger distance between test and each row of the csr sparse matrix

有什么方法可以使用 numpy 符号返回所需的距离向量，也许使用 np.apply_along_axis方法？我以前见过这样做，但似乎无法在这里得到它。提前致谢。

注意:我想避免显式的 for 循环，因为这样效率很低。我正在寻找最优化/最快的方法来做到这一点。

Answer 1

矢量化解决方案

这是我通过一些优化和一个关键技巧得出的最终矢量化解决方案，假设 s 作为 csr_matrix 类型的输入稀疏矩阵。

k1 = np.sqrt(1/2)
k2s = np.sqrt(test.dot(test))
out = k1*np.sqrt(k2s + s.sum(1).A1 -2*np.sqrt(s*test))

---

回放历史

经过一系列优化后，最终的矢量化解决方案已经达到，我会尝试回放以供我和其他人引用，对于在这里冗长，我深表歉意，但我觉得这是必要的。

第 1 阶段

从循环中的函数定义开始:

N = s.shape[0]
out = np.zeros(N)
for i in range(s.shape[0]):
    ai = s[i].toarray()
    out[i] = np.sqrt(1/2) * np.sqrt(  np.sum((np.sqrt(test) - np.sqrt(ai))**2)  )

第 2 阶段

取出常量并在外面计算平方根:

k1 = np.sqrt(1/2)
k2 - np.sqrt(test)

N = s.shape[0]
out = np.zeros(N)
for i in range(s.shape[0]):
    ai = s[i].toarray()
    out[i] = np.sum((k2 - np.sqrt(ai))**2)

out = np.sqrt(out)
out *= k1

阶段 #3(关键技巧)

这里的关键技巧是我们将使用数学公式:

(A-B)**2 = A**2) + B**2 - 2*A*B

因此，

sum((A-B)**2) = sum(A**2) + sum(B**2) - 2*sum(A*B)

最后一部分 sum(A*B) 只是矩阵乘法，这是这里的主要性能提升器。

简化为:

k1 = np.sqrt(1/2)
k2 - np.sqrt(test)

N = s.shape[0]
out = np.zeros(N)
for i in range(s.shape[0]):
    ai = s[i].toarray()
    out[i] = (k2**2).sum() + (np.sqrt(ai))**2).sum() -2*np.sqrt(ai).dot(k2)

out = np.sqrt(out)
out *= k1

进一步简化为:

k1 = np.sqrt(1/2)
k2 - np.sqrt(test)

N = s.shape[0]
out = np.zeros(N)
for i in range(s.shape[0]):
    ai = s[i].toarray()
    out[i] = (k2**2).sum() + ai.sum() -2*np.sqrt(ai).dot(k2)

out = np.sqrt(out)
out *= k1

第 4 阶段

获取常量 (k2**2).sum() 并获取稀疏矩阵的逐行求和:

k1 = np.sqrt(1/2)
k2 - np.sqrt(test)
k2s = (k2**2).sum()

N = s.shape[0]
out = np.zeros(N)
for i in range(s.shape[0]):
    ai = s[i].toarray()
    out[i] =  -2*np.sqrt(ai).dot(k2)

out += k2s + s.sum(1).A1 # row-wise summation of sparse matrix added here
out = np.sqrt(out)
out *= k1

第 5 阶段

最后的技巧是完全删除循环。因此，在循环中，我们使用 np.sqrt(s[i]).dot(k2) 计算每个输出元素。矩阵乘法可以简单地在所有行上完成:np.sqrt(s)*k2。就这样!

剩下的将是:

k1 = np.sqrt(1/2)
k2 - np.sqrt(test)
k2s = (k2**2).sum()

out = -2*np.sqrt(s)*k2 # Loop gone here
out += k2s + s.sum(1).A1
out = np.sqrt(out)
out *= k1

在使用 inner 点积得到 k2s 之后简化为 -

k1 = np.sqrt(1/2)
k2 = np.sqrt(test)
k2s = k2.dot(k2)
out = k1*np.sqrt(k2s + s.sum(1).A1 -2*np.sqrt(s)*k2)

我们可以避免 test 的平方根计算来得到 k2 从而进一步简化像这样的事情 -

k1 = np.sqrt(1/2)
k2s = np.sqrt(test.dot(test))
out = k1*np.sqrt(k2s + s.sum(1).A1 -2*np.sqrt(s*test))