Jar*_*edL 3 python numpy r normal-distribution matrix
编辑2:这篇文章似乎已经从CrossValidated转移到StackOverflow,因为它主要是关于编程,但这意味着花哨的MathJax不再起作用了.希望这仍然可读.
假设我想计算两个向量之间的平方Mahalanobis距离x和y协方差矩阵S.这是一个相当简单的函数定义
M2(x, y; S) = (x - y)^T * S^-1 * (x - y)
使用python的numpy包我可以这样做
# x, y = numpy.ndarray of shape (n,)
# s_inv = numpy.ndarray of shape (n, n)
diff = x - y
d2 = diff.T.dot(s_inv).dot(diff)
或者在R中
diff <- x - y
d2 <- t(diff) %*% s_inv %*% diff
但就我而言,我得到了
m按n矩阵Xn- 三维矢量 mun通过n协方差矩阵S并希望找到m维向量d,使得
d_i = M2(x_i, mu; S) ( i = 1 .. m )
哪排x_i是第一i排X.
使用python中的简单循环来实现这一点并不困难:
d = numpy.zeros((m,))
for i in range(m):
diff = x[i,:] - mu
d[i] = diff.T.dot(s_inv).dot(diff)
当然,假设外部循环发生在python而不是numpy库中的本机代码意味着它没有尽可能快.$ n $和$ m $分别约为3-4和几十万,我在交互式程序中经常这样做,所以加速会非常有用.
在数学上,我能够使用基本矩阵运算来表达这一点的唯一方法是
d = diag( X' * S^-1 * X'^T )
哪里
x'_i = x_i - mu
编写矢量化版本很简单,但遗憾的是计算一个100亿以上的元素矩阵并且只采用对角线的效率低得多...我相信这个操作应该可以使用爱因斯坦符号很容易表达,因此希望能与快速评估numpy的einsum功能,但我还没有开始弄清楚,黑魔法是如何工作的.
所以,我想知道:有没有更好的方法以数学方式(在简单的矩阵运算方面)制定这个操作,或者有人可以提出一些很好的矢量化(python或R)代码来有效地做到这一点？
我实际上并不想这样做一次,我想这样做k~100次.鉴于:
m按n矩阵X
k按n矩阵U
设置的n由n协方差矩阵的每个表示S_j(j = 1..k)
m通过k矩阵找到D这样的
D_i,j = M(x_i, u_j; S_j)
其中i = 1..m,j = 1..k,x_i是i的第i行X和u_j是j的第i行U.
即,矢量化以下代码:
# s_inv is (k x n x n) array containing "stacked" inverses
# of covariance matrices
d = numpy.zeros( (m, k) )
for j in range(k):
for i in range(m):
diff = x[i, :] - u[j, :]
d[i, j] = diff.T.dot(s_inv[j, :, :]).dot(diff)
首先,似乎你可能正在获得S然后反转它.你不应该这样做; 它的速度慢,数值不准确.相反,你应该得到S的Cholesky因子L,使得S = LL ^ T; 然后
M^2(x, y; L L^T)
= (x - y)^T (L L^T)^-1 (x - y)
= (x - y)^T L^-T L^-1 (x - y)
= || L^-1 (x - y) ||^2,
由于L是三角形,因此可以有效地计算L ^ -1(x-y).
事实证明,scipy.linalg.solve_triangular如果你正确地重塑它,我们会很高兴地同时做一些这些:
L = np.linalg.cholesky(S)
y = scipy.linalg.solve_triangular(L, (X - mu[np.newaxis]).T, lower=True)
d = np.einsum('ij,ij->j', y, y)
稍微打破一下,y[i, j]是L ^ -1(X_j - \mu)的第i个分量.然后einsum打电话
d_j = \sum_i y_{ij} y_{ij}
= \sum_i y_{ij}^2
= || y_j ||^2,
我们需要的.
不幸的是,solve_triangular它不会在第一个参数上进行矢量化,所以你应该只是在那里循环.如果k只有大约100,那就不会是一个重大问题.
如果你实际上给了S ^ -1而不是S,那么你确实可以einsum更直接地做到这一点.由于S在您的情况下非常小,因此实际上反转矩阵也可能会更快.但是,只要n是一个非常重要的大小,你就会通过这样做而丢掉很多数值精度.
要弄清楚如何处理einsum,请根据组件编写所有内容.我将直接进入奖金案例,为了方便起见,写下S_j ^ -1 = T_j:
D_{ij} = M^2(x_i, u_j; S_j)
= (x_i - u_j)^T T_j (x_i - u_j)
= \sum_k (x_i - u_j)_k ( T_j (x_i - u_j) )_k
= \sum_k (x_i - u_j)_k \sum_l (T_j)_{k l} (x_i - u_j)_l
= \sum_{k l} (X_{i k} - U_{j k}) (T_j)_{k l} (X_{i l} - U_{j l})
因此,如果我们制作X形状(m, n),U形状(k, n)和T形状的数组(k, n, n),那么我们可以将其写成
diff = X[np.newaxis, :, :] - U[:, np.newaxis, :]
D = np.einsum('jik,jkl,jil->ij', diff, T, diff)
哪里diff[j, i, k] = X_[i, k] - U[j, k].