如何用Numpy实现向量元素乘积矩阵以计算误差函数V？

问题描述 #

我需要计算误差函数V：

V = Σi Σj X[i] X[j] σ[i][j]

其中σ[i][j]是给定矩阵，希望找到高效的Numpy实现方案。原本计划创建矩阵Y[i][j] = X[i]*X[j]，再对Y * σ求和，但之前尝试用meshgrid(X, X)后对每行应用np.prod，不仅没得到预期结果，还需要Python循环，效率低下。

补充可复现代码示例：

import numpy as np

cov = np.array(((0.1, 0.05), (0.05, 0.25)))
x = np.array((0.6, 0.4))

# 预期计算结果
desired = x[0]*x[0]*cov[0][0] + x[0]*x[1]*cov[0][1] + x[1]*x[0]*cov[1][0] + x[1]*x[1] * cov[1][1]

解决方案 #

最优方案：直接利用二次型矩阵运算 #

你要计算的V本质是二次型，可以通过Numpy的矩阵乘法一步完成，完全不需要额外构造Y矩阵，内存和计算效率最优：

v = x.T @ cov @ x
print(v)  # 输出：0.1，与desired结果一致

备选方案：构造Y矩阵后求和 #

如果坚持要生成Y[i][j] = X[i]*X[j]，可以用Numpy的外积函数np.outer快速生成Y，再逐元素相乘求和：

y = np.outer(x, x)
v = np.sum(y * cov)
print(v)  # 输出：0.1

修正meshgrid的用法 #

如果想用meshgrid实现，不需要循环，直接对生成的两个网格数组逐元素相乘即可得到Y矩阵：

x1, x2 = np.meshgrid(x, x)
y = x1 * x2
v = np.sum(y * cov)
print(v)  # 输出：0.1

之前的错误在于用np.prod处理每行，而meshgrid返回的两个数组已经对应了所有X[i]和X[j]的组合，直接相乘就能得到正确的Y矩阵。

内容的提问来源于stack exchange，提问作者David K.