努比半正定预警

Question

在我编写的一个Python脚本中，我用表达式模拟多变量正态随机向量

np.random.multivariate_normal(np.zeros(dim_obs), y_cov)

我的脚本运行，但生成以下警告：

RuntimeWarning: covariance is not positive-semidefinite.

还有我在这里输入的一些调试打印语句，大多数时候都会打印False

print( np.all(np.linalg.eigvals(y_cov) > 0) )

为什么会有误报？我的y_cov是半正定的，因为它是（很抱歉没有TeX标记）B x x'B'+y'，其中B是一个矩阵，其他的是每个元素都是正的随机向量。

在这个特殊的运行B实际上只是一个大小为9的1向量。我能忽略这个警告吗？From the documentation：

Note that the covariance matrix must be positive semidefinite (a.k.a. nonnegative-definite). Otherwise, the behavior of this method is undefined and backwards compatibility is not guaranteed.

编辑： 这是一个完全可以运行的东西。谢谢你的提示@user2357112。

import numpy as np
num_factors = 1
dim_obs = 9
u = np.random.normal(size = num_factors)
v = np.random.normal(size = dim_obs)
y_cov = np.dot(np.ones((9,1)), np.exp(u.reshape((num_factors,1))/2))
y_cov = np.dot(y_cov, np.exp(u.reshape((1,num_factors))/2)) #transpose
y_cov = np.dot(y_cov, np.transpose(np.ones((9,1))))
y_cov += np.dot(np.exp( v.reshape((dim_obs,1)) / 2), 
                np.exp( v.reshape((1,dim_obs)) / 2))
print( np.random.multivariate_normal(np.zeros(dim_obs), y_cov) )
print( np.all(np.linalg.eigvals(y_cov) > 0) )
print( np.linalg.eigvals(y_cov)  )

Answer 1

理论上，你的矩阵是半正定的，几个特征值正好为零。但浮点数的计算引入了截断误差，导致某些特征值很小，但却负；因此，矩阵不是半正定的。

目前看来，警告可能会被忽略；但是NumPy文档说，在非psd情况下的行为是未定义的，所以我不想依赖于此。纠正浮点错误的一种方法是在y_cov中添加一个很小的单位矩阵倍数。例如，如下所示：

min_eig = np.min(np.real(np.linalg.eigvals(y_cov)))
if min_eig < 0:
    y_cov -= 10*min_eig * np.eye(*y_cov.shape)

增加一个固定的恒等式倍数，如1e-12，对所有合理的大小矩阵都有效，但对结果仍然无关紧要。

---

为了完整起见，复制问题的一种更简单的方法是：

import numpy as np
x = np.random.normal(size=(5,))
y = np.outer(x, x)
z = np.random.multivariate_normal(np.zeros(5), y)    

这会抛出同样的警告（概率很高）。

Answer 2

在您的例子中，生成高斯样本的一种更有效的方法是观察协方差矩阵等于a*a.T + b*b.T（a）的多变量、零均值高斯随机向量，它也不受@zaq所识别的数值问题的影响，b：列向量）在分布上与随机向量a*w1 + b*w2相等，其中w1和w2是零均值和方差的独立高斯标量随机变量1