在Python中计算log(det(AA^T)+1)时出现数学错误。

Question

我正在尝试在Python中估算log(det(AA^T)+1)的平均值。我的简单代码在处理17×17矩阵时出现了数学错误。以下是代码：

iter = 10000
for n in xrange(1,20):
    h = n
    dets = []
    for _ in xrange(iter):
        A = (np.random.randint(2, size=(h,n)))*2-1
        detA_Atranspose = np.linalg.det(np.dot(A, A.transpose()))
        try:
            logdetA_Atranspose = math.log(detA_Atranspose+1,2)
        except ValueError:
            print "Ooops!", n,detA_Atranspose
        dets.append(logdetA_Atranspose)
    print np.mean(dets)

A应该是一个元素为-1或1的矩阵。

我做错了什么，怎么修复？17有什么特别之处吗？

- Simd

---

17没有什么特别的地方，只是因为你无法计算log detA_Atranspose <= 0。 - jtitusj

@JohnTitusJungao 为什么在n < 17时这些矩阵不存在？ - Simd

我在 n = 16 时遇到了错误。这里是截图。 - jtitusj

@JohnTitusJungao 哦 :) - Simd

你应该将标题更正为实际使用的公式 log(1+det(AA^T), 2)。 - Lutz Lehmann

Answer 1

对于标题中的公式（之前是logdet(AA^T)）：

对于某些随机矩阵A，det(AA^T) 可能简单地等于0。因为计算 log(0) 是无效的，所以函数会失败。

请注意，在理论上，det(AA^T) 不能为负，因为AA^T是一个半正定矩阵（这意味着所有特征值都是非负的，并且意味着det>=0）。

对于代码中的公式（logdet(1+AA^T)）

您应该使用numpy.linalg.slogdet()并计算slogdet(1+A.dot(A.T))

从其文档中得知：

"计算数组的符号和（自然）对数行列式。

如果一个数组的行列式非常小或非常大，那么调用det可能会溢出或下溢。这个程序更加健壮，因为它计算的是行列式的对数，而不是行列式本身。

- Tomer Levinboim

---

你对 log(1+AA^T) 中的 +1 是正确的。我想我是从标题中拿来的...所以，这只能是由于溢出而发生的。 - Tomer Levinboim

建议使用slogdet应该将计算转换为浮点数，这也将修复整数溢出问题。然而，有相当小的概率，行列式仍可能为零。 - Lutz Lehmann

"det(AA^T) 不能为负数" ...但根据代码似乎是负数。 - Simd

这可能是您原始问题的另一种解决方案：计算logdet(1 + AA^T + pI)，其中p是一个小标量（比如1e-6），I是单位矩阵。这应该可以消除不稳定性问题，但会牺牲一些准确性。 - Tomer Levinboim

@LutzL 你确定吗？因为在64位有符号整数上，溢出会发生在2^63处，而 ln(2^63) = 43.6，已经观察到的值并不是不可能的...我不确定 numpy.linalg.det 的内部是使用浮点数还是整数运算，但文档确实说它容易发生溢出。 - Cimbali