python - 乔列斯基分解浮点误差

Question

在不使用 numpy.linalg.* 的情况下在 Python 中实现 Choleky 分解(是的，这是一个作业)我遇到了有关使用 float 的问题。我的算法适用于常规整数:

    for i in range(n):
        for k in range(i + 1):
            sum = 0
            for j in range(k):
                sum += L[i, j] * L[k, j]
            if i == k:
                if M[i, i] - sum < 0:
                    raise ValueError("Matrix is non-positive definite")
                L[i, i] = np.sqrt(M[i, i] - sum)
            else:
                if np.isclose(L[k, k] * (M[i, k] - sum), 0):
                    raise ValueError("Matrix is non-positive definite")
                L[i, k] = (L[k, k] * (M[i, k] - sum)) ** (-1)

我事先测试了矩阵的对称性； n 是维度，L 成为下三角 Cholesky 因子。

使用随机浮点 nxn 矩阵乘以它们的转置(为了获得正定矩阵)，两个 ValueErrors 都会被引发，即，如果不提高 ValueErrors，L 输出矩阵将部分填充 NaN 和 inf 值。如何在 python 中使用 float ？

编辑:最小可重复示例:

M = np.random.randn(2, 2)
M = np.dot(M, M.transpose())
# produces for example [[0.68283219, -0.33497034], [-0.33497034, 0.72113541]]
run_cholesky(M)

Answer 1

将 M[i, k] 保存在变量中，然后进行减法而不是求和即可解决问题:

for i in range(n):
    for k in range(i + 1):
        val = M[i, k]
        for j in range(k):
            val -= L[i, j] * L[k, j]
        if i == k:
            if val < 0:
                raise ValueError("Matrix is non-positive definite")
            L[i, k] = np.sqrt(val)
        else:
            if np.isclose(L[k, k], 0):
                raise ValueError("Matrix is non-positive definite")
            L[i, k] = val / L[k, k]