我觉得你就快成功了。你有一个矩阵A = np.squeeze((X[0] - mu)[:, None] @ S)，你把它和B = ((X[0] - mu)).T)相乘，但是你只需要对角元素。
正如所指出的，hereC = 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)

这也可以用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)))