cholesky分解浮点错误

Question

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

    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 都被引发，即 w/out 引发 ValueErrors L 输出矩阵部分填充 NaN 和 inf值。如何在 python 中使用浮点数？

编辑：最小可重现示例：

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 = np.zeros((n, n), dtype=float)

Answer 2

将 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]

cholesky分解浮点错误

cholesky decomposition floating point error

python

floating-point

numpy

matrix-decomposition