英文:
Broadcasting for numpy array - vectorized quadratic form
问题
I'd like to compute the part of multivariate normal distribution density that is a quadratic form:
要计算多元正态分布密度的二次形式部分：
(X - mu)^T * S * (X - mu)
Assume the data:
假设数据如下：
mu = np.array([[1,2,3], [4,5,6]])
S = np.array([np.eye(3)*3, np.eye(3)*5])
X = np.array([np.random.random(3*10)]).reshape(10, 3)
Now, an iterative process would be to calculate:
现在，一个迭代的方法是计算：
(X[0] - mu[0]) @ S[0] @ (X[0] - mu[0]).T, (X[0] - mu[1]) @ S[1] @ (X[0] - mu[1]).T
However, I guess that's not the fastest approach. What I tried is:
不过，我猜这不是最快的方法。我尝试的是：
np.squeeze((X[0] - mu)[:, None] @ S) @ ((X[0] - mu)).T
But the values that I want are placed on the main diagonal of the matrix above. I could use np.diagonal(), but is there a better way to perform the calculations?
但我想要的值位于上面矩阵的主对角线上。我可以使用 np.diagonal()，但是否有更好的方法来执行计算？
英文:
I'd like to compute the part of multivariate normal distribution density that is a quadratic form
(X - mu)^T * S * (X - mu)
Assume the data
mu = np.array([[1,2,3], [4,5,6]])
S = np.array([np.eye(3)*3, np.eye(3)*5])
X = np.array([np.random.random(3*10)]).reshape(10, 3)
Now, an iterative process would be to calculate
(X[0] - mu[0]) @ S[0] @ (X[0] - mu[0]).T, (X[0] - mu[1]) @ S[1] @ (X[0] - mu[1]).T
(I don't need to vectorize with respect to X). However, I guess that's not the fastest approach. What I tried is
np.squeeze((X[0] - mu)[:, None] @ S) @ ((X[0] - mu)).T
But the values that I want are placed on the main diagonal of matrix above. I could use np.diagonal(), but is there a better way to perform the calculations?
答案1
得分: 2
我认为你差不多已经做到了。你有矩阵 A = np.squeeze((X[0] - mu)[:, None] @ S)，然后你对它进行了矩阵乘法 B = ((X[0] - mu)).T)，但你只想要对角元素。
正如在这里指出的那样，C = N.diag(A.dot(B)) 等同于 C = (A * B.T).sum(-1)，这导致了以下解决方案：
import numpy as np
mu = np.array([[1,2,3], [4,5,6]])
S = np.array([np.eye(3)*3, np.eye(3)*5])
X = np.array([np.random.random(3*10)]).reshape(10, 3)
res1 = (X[0] - mu[0]) @ S[0] @ (X[0] - mu[0]).T, (X[0] - mu[1]) @ S[1] @ (X[0] - mu[1]).T
res2 = (np.squeeze((X[0] - mu)[:, None] @ S) * (X[0] - mu)).sum(-1)
print(res1)
print(res2)
英文:
I think you were almost there. You have matrix A = np.squeeze((X[0] - mu)[:, None] @ S) which you matrix multiply with B = ((X[0] - mu)).T) but you only want the diagonal elements.
As pointed out here C = N.diag(A.dot(B)) is equivalent to C = (A * B.T).sum(-1) which leads to the following solution:
import numpy as np
mu = np.array([[1,2,3], [4,5,6]])
S = np.array([np.eye(3)*3, np.eye(3)*5])
X = np.array([np.random.random(3*10)]).reshape(10, 3)
res1 = (X[0] - mu[0]) @ S[0] @ (X[0] - mu[0]).T, (X[0] - mu[1]) @ S[1] @ (X[0] - mu[1]).T
res2 = (np.squeeze((X[0] - mu)[:, None] @ S) * (X[0] - mu)).sum(-1)
print(res1)
print(res2)
答案2
得分: 1
这也可以使用 np.einsum 表达，允许您在 X 上进行广播：
import numpy as np
mu = np.array([[1,2,3], [4,5,6]])
S = np.array([np.eye(3)*3, np.eye(3)*5])
X = np.random.random((10, 3))
resOP = np.array([(X[0] - mu[0]) @ S[0] @ (X[0] - mu[0]).T, (X[0] - mu[1]) @ S[1] @ (X[0] - mu[1]).T])
resNin17 = np.einsum("...ij, ...ij->...i", np.einsum("...j, ...ij", (X[:, None] - mu), S), (X[:, None] - mu))
assert np.allclose(resOP, resNin17[0])
或者只计算一行：
assert np.array_equal(resNin17[2], np.einsum("...ij, ...ij->...i", np.einsum("...j, ...ij", (X[2] - mu), S), (X[2] - mu)))
英文:
This can also be expressed using np.einsum allowing you to broadcast over X as well:
import numpy as np
mu = np.array([[1,2,3], [4,5,6]])
S = np.array([np.eye(3)*3, np.eye(3)*5])
X = np.random.random((10, 3))
resOP = np.array([(X[0] - mu[0]) @ S[0] @ (X[0] - mu[0]).T, (X[0] - mu[1]) @ S[1] @ (X[0] - mu[1]).T])
resNin17 = np.einsum("...ij, ...ij->...i", np.einsum("...j, ...ij", (X[:, None] - mu), S), (X[:, None] - mu))
assert np.allclose(resOP, resNin17[0])
Or to just calculate one row:
assert np.array_equal(resNin17[2], np.einsum("...ij, ...ij->...i", np.einsum("...j, ...ij", (X[2] - mu), S), (X[2] - mu)))
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论