确保正定协方差矩阵答案

【问题标题】：Ensuring positive definite covariance matrix确保正定协方差矩阵
【发布时间】：2017-03-27 04:32:26
【问题描述】：

我的神经网络的输出充当协方差矩阵的条目。但是，输出和条目之间的一一对应导致非正定协方差矩阵。

因此，我读到了https://www.quora.com/When-carrying-out-the-EM-algorithm-how-do-I-ensure-that-the-covariance-matrix-is-positive-definite-at-all-times-avoiding-rounding-issues 和https://en.wikipedia.org/wiki/Cholesky_decomposition，更具体地说，“当 A 有实条目时，L 也有实条目，分解可以写成 A = LL^T”。

现在我的输出对应于 L 矩阵的条目，然后我通过将其乘以它的转置来生成协方差矩阵。

但是，有时我仍然会遇到非正定矩阵的错误。这怎么可能？

我发现了一个产生错误的矩阵，请参阅

print L.shape
print Sigma.shape

S = Sigma[1,18,:,:] # The matrix that gives the error
L_ = L[1,18,:,:]
print L_
S = np.dot(L_,np.transpose(L_))
print S
chol = np.linalg.cholesky(S)

作为输出：

(3, 20, 2, 2)
(3, 20, 2, 2)
[[ -1.69684255e+00   0.00000000e+00]
 [ -1.50235415e+00   1.73807144e-04]]
[[ 2.87927461  2.54925847]
 [ 2.54925847  2.25706792]]
.....
LinAlgError: Matrix is not positive definite

但是，复制值的这段代码可以正常工作（但可能不是完全相同的值，因为并非所有小数都被打印）

B = np.array([[-1.69684255e+00, 0.00000000e+00], [-1.50235415e+00, 1.73807144e-04]])
A = np.dot(B,B.T)
chol_A = np.linalg.cholesky(A)

所以问题是：

使用 Sigma = LL' 的方法是否正确（使用 ' 转置）？
如果是，为什么会出现错误？这可能是由于舍入问题造成的吗？

编辑：我还计算了特征值

print np.linalg.eigvalsh(S)
[ -7.89378944432428397703915834426880e-08
   5.13634252548217773437500000000000e+00]

对于第二种情况

print np.linalg.eigvalsh(A)
[  1.69341869415973178547574207186699e-08
   5.13634263409323210680668125860393e+00]

所以第一种情况有一个轻微的负特征值，它声明了非正定性。但是如何解决呢？

【问题讨论】：

好吧，您引用的Cholesky decomposition 仅针对作为PD 的矩阵S 定义。句子“当 A 有真实条目时，L 也有真实条目，并且分解可以写成 A = LL^T”假设 A 是 PD，而您的 S 显然不是，正如您所观察到的。跨度>
如果您将 L 的值打印得更精确并复制它们，您是否会在原件和副本之间得到相同的差异？
是的，np.set_printoptions(precision=40) 的结果相同。通过向协方差矩阵添加 eps*I 来解决它，尽管它似乎不是最好的解决方案
你能解决这个问题吗？我正在努力实现同样的目标

标签： numpy neural-network tensorflow deep-learning linear-algebra

【解决方案1】：

这看起来像是一个数值问题，但一般来说，LL' 永远是正定的（如果 L 是可逆的）是不正确的。例如，将 L 作为矩阵，其中每列为 [1 0 0 0 ... 0]（或更极端 - 将 L 视为任意维数的零矩阵），LL' 不会是 PD。一般来说，我会建议这样做

S = LL' + eps I

它同时解决了这两个问题（对于小 eps），并且是一个“正则化”协方差估计。您甚至可以使用 Ledoit-Wolf 估计器获得 eps 的“最佳”（在某些假设下）值。

【讨论】：

【解决方案2】：

我怀疑L*L' 的计算在第一种情况下是用浮点数完成的，在第二种情况下是用双精度数完成的。我尝试将你的 L 作为浮点矩阵，计算 L*L' 并找到它的特征值，我得到的值与你在第一种情况下所做的相同，但如果我将 L 转换为双精度矩阵，计算 L*L' 和找到与第二种情况相同的特征值。

这是有道理的，因为在计算 L*L'[1,1] 时，浮点数中 1.73807144e-04 的平方与 -1.50235415e+00 的平方相比可以忽略不计。

如果我是对的，解决方案是在任何计算之前将 L 转换为双精度矩阵。

【讨论】：