用于计算（平方）马氏距离的矢量化代码
Vectorizing code to calculate (squared) Mahalanobis Distiance
编辑 2：此 post 似乎已从 CrossValidated 移至 Whosebug，因为它主要与编程有关，但这意味着花哨的 MathJax 不再起作用。希望这仍然是可读的。
假设我想用协方差矩阵 S 计算两个向量 x 和 y 之间的平方马氏距离。这是一个由
定义的相当简单的函数
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 by n 协方差矩阵 S
并且想要找到m维向量d使得
d_i = M2(x_i, mu; S) ( i = 1 .. m )
其中 x_i 是 X 的第 i 行。
使用 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)
通过 k 矩阵 D 找到 m 使得
D_i,j = M(x_i, u_j; S_j)
其中 i = 1..m、j = 1..k、x_i 是 X 的第 i 行，而 u_j 是 [=59=第 ] 行 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 = L L^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})
因此，如果我们制作形状为 (m, n) 的数组 X、形状为 (k, n) 的数组 U 以及形状为 (k, n, n) 的数组 T , 那么我们可以写成
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].
Dougal 以出色而详细的答案解决了这个问题，但我认为我会分享一个小的修改，我发现它可以提高效率，以防其他人试图实现它。开门见山：
Dougal 的方法如下：
def mahalanobis2(X, mu, sigma):
L = np.linalg.cholesky(sigma)
y = scipy.linalg.solve_triangular(L, (X - mu[np.newaxis,:]).T, lower=True)
return np.einsum('ij,ij->j', y, y)
我试过的一个数学上等效的变体是
def mahalanobis2_2(X, mu, sigma):
# Cholesky decomposition of inverse of covariance matrix
# (Doing this in either order should be equivalent)
linv = np.linalg.cholesky(np.linalg.inv(sigma))
# Just do regular matrix multiplication with this matrix
y = (X - mu[np.newaxis,:]).dot(linv)
# Same as above, but note different index at end because the matrix
# y is transposed here compared to above
return np.einsum('ij,ij->i', y, y)
运行 两个版本使用相同的 运行dom 输入并记录时间（以毫秒为单位）进行 20 倍的正面交锋。对于 X 作为 1,000,000 x 3 矩阵（mu 和 sigma 3 和 3x3），我得到：
Method 1 (min/max/avg): 30/62/49
Method 2 (min/max/avg): 30/47/37
第二个版本的速度大约提高了 30%。我主要是 运行 在 3 或 4 维中这样做，但要看看它是如何缩放的我尝试将 X 设为 1,000,000 x 100 并得到：
Method 1 (min/max/avg): 970/1134/1043
Method 2 (min/max/avg): 776/907/837
大致相同的改进。
我在对 Dougal 的回答的评论中提到了这一点，但在此处添加以提高知名度：
上面的第一对方法采用单个中心点 mu 和协方差矩阵 sigma 并计算到 X 的每一行的平方马氏距离。我的奖励问题是多次执行此操作有很多组 mu 和 sigma 并输出一个二维矩阵。上面的一组方法可用于通过简单的 for 循环来完成此操作，但 Dougal 还使用 einsum.
发布了一个更聪明的示例
我决定通过使用它们解决以下问题来比较这些方法：给定 k d 维正态分布（中心存储在 k 行中通过d矩阵U和kbydbyd数组S最后两个维度的协方差矩阵，求密度n 点存储在 n by d 矩阵 X.
的行中
多元正态分布的密度是点与均值的平方马氏距离的函数。 Scipy 有一个作为 scipy.stats.multivariate_normal.pdf 的实现以供参考。我 运行 所有三种方法每次都使用相同的 运行dom 参数相互对抗 10 倍，d=3, k=96, n=5e5。以下是结果，在 points/sec:
[Method]: (min/max/avg)
Scipy: 1.18e5/1.29e5/1.22e5
Fancy 1: 1.41e5/1.53e5/1.48e5
Fancy 2: 8.69e4/9.73e4/9.03e4
Fancy 2 (cheating version): 8.61e4/9.88e4/9.04e4
其中 Fancy 1 是上述两种方法中较好的一种，Fancy2 是 Dougal 的第二种解决方案。由于 Fancy 2 需要计算所有协方差矩阵的逆，我还尝试了 "cheating version" 并将其作为参数传递，但看起来并没有什么不同。我曾计划包括非矢量化实施，但那太慢了，可能要花一整天的时间。
我们可以从中得出的结论是，使用 Dougal 的第一种方法比 Scipy 快 20%。 不幸的是，尽管第二种方法很聪明只比第一个快 60%。可能还有一些其他的优化可以完成，但这对我来说已经足够快了。
我还测试了它如何随着更高的维度进行缩放。随着 d=100, k=96, n=1e4:
Scipy: 7.81e3/7.91e3/7.86e3
Fancy 1: 1.03e4/1.15e4/1.08e4
Fancy 2: 3.75e3/4.10e3/3.95e3
Fancy 2 (cheating version): 3.58e3/4.09e3/3.85e3
Fancy 1这次似乎优势更大。同样值得注意的是 Scipy 抛出了 LinAlgError 8/10 次，可能是因为我的一些 运行domly 生成的 100x100 协方差矩阵接近奇异（这可能意味着其他两种方法数值上不太稳定，我没有实际检查结果）。
编辑 2：此 post 似乎已从 CrossValidated 移至 Whosebug，因为它主要与编程有关，但这意味着花哨的 MathJax 不再起作用。希望这仍然是可读的。
假设我想用协方差矩阵 S 计算两个向量 x 和 y 之间的平方马氏距离。这是一个由
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维向量munbyn协方差矩阵S
并且想要找到m维向量d使得
d_i = M2(x_i, mu; S) ( i = 1 .. m )
其中 x_i 是 X 的第 i 行。
使用 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矩阵Xk通过n矩阵U由
n个协方差矩阵组成的n集，每个表示为S_j(j = 1..k)
通过 k 矩阵 D 找到 m 使得
D_i,j = M(x_i, u_j; S_j)
其中 i = 1..m、j = 1..k、x_i 是 X 的第 i 行，而 u_j 是 [=59=第 ] 行 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 = L L^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})
因此，如果我们制作形状为 (m, n) 的数组 X、形状为 (k, n) 的数组 U 以及形状为 (k, n, n) 的数组 T , 那么我们可以写成
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].
Dougal 以出色而详细的答案解决了这个问题，但我认为我会分享一个小的修改，我发现它可以提高效率，以防其他人试图实现它。开门见山：
Dougal 的方法如下：
def mahalanobis2(X, mu, sigma):
L = np.linalg.cholesky(sigma)
y = scipy.linalg.solve_triangular(L, (X - mu[np.newaxis,:]).T, lower=True)
return np.einsum('ij,ij->j', y, y)
我试过的一个数学上等效的变体是
def mahalanobis2_2(X, mu, sigma):
# Cholesky decomposition of inverse of covariance matrix
# (Doing this in either order should be equivalent)
linv = np.linalg.cholesky(np.linalg.inv(sigma))
# Just do regular matrix multiplication with this matrix
y = (X - mu[np.newaxis,:]).dot(linv)
# Same as above, but note different index at end because the matrix
# y is transposed here compared to above
return np.einsum('ij,ij->i', y, y)
运行 两个版本使用相同的 运行dom 输入并记录时间（以毫秒为单位）进行 20 倍的正面交锋。对于 X 作为 1,000,000 x 3 矩阵（mu 和 sigma 3 和 3x3），我得到：
Method 1 (min/max/avg): 30/62/49
Method 2 (min/max/avg): 30/47/37
第二个版本的速度大约提高了 30%。我主要是 运行 在 3 或 4 维中这样做，但要看看它是如何缩放的我尝试将 X 设为 1,000,000 x 100 并得到：
Method 1 (min/max/avg): 970/1134/1043
Method 2 (min/max/avg): 776/907/837
大致相同的改进。
我在对 Dougal 的回答的评论中提到了这一点，但在此处添加以提高知名度：
上面的第一对方法采用单个中心点 mu 和协方差矩阵 sigma 并计算到 X 的每一行的平方马氏距离。我的奖励问题是多次执行此操作有很多组 mu 和 sigma 并输出一个二维矩阵。上面的一组方法可用于通过简单的 for 循环来完成此操作，但 Dougal 还使用 einsum.
我决定通过使用它们解决以下问题来比较这些方法：给定 k d 维正态分布（中心存储在 k 行中通过d矩阵U和kbydbyd数组S最后两个维度的协方差矩阵，求密度n 点存储在 n by d 矩阵 X.
多元正态分布的密度是点与均值的平方马氏距离的函数。 Scipy 有一个作为 scipy.stats.multivariate_normal.pdf 的实现以供参考。我 运行 所有三种方法每次都使用相同的 运行dom 参数相互对抗 10 倍，d=3, k=96, n=5e5。以下是结果，在 points/sec:
[Method]: (min/max/avg)
Scipy: 1.18e5/1.29e5/1.22e5
Fancy 1: 1.41e5/1.53e5/1.48e5
Fancy 2: 8.69e4/9.73e4/9.03e4
Fancy 2 (cheating version): 8.61e4/9.88e4/9.04e4
其中 Fancy 1 是上述两种方法中较好的一种，Fancy2 是 Dougal 的第二种解决方案。由于 Fancy 2 需要计算所有协方差矩阵的逆，我还尝试了 "cheating version" 并将其作为参数传递，但看起来并没有什么不同。我曾计划包括非矢量化实施，但那太慢了，可能要花一整天的时间。
我们可以从中得出的结论是，使用 Dougal 的第一种方法比 Scipy 快 20%。 不幸的是，尽管第二种方法很聪明只比第一个快 60%。可能还有一些其他的优化可以完成，但这对我来说已经足够快了。
我还测试了它如何随着更高的维度进行缩放。随着 d=100, k=96, n=1e4:
Scipy: 7.81e3/7.91e3/7.86e3
Fancy 1: 1.03e4/1.15e4/1.08e4
Fancy 2: 3.75e3/4.10e3/3.95e3
Fancy 2 (cheating version): 3.58e3/4.09e3/3.85e3
Fancy 1这次似乎优势更大。同样值得注意的是 Scipy 抛出了 LinAlgError 8/10 次，可能是因为我的一些 运行domly 生成的 100x100 协方差矩阵接近奇异（这可能意味着其他两种方法数值上不太稳定，我没有实际检查结果）。