3

在不使用 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 都会被提升，即在没有提升 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)
4

2 回答 2

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]
于 2019-12-15T12:47:56.533 回答
0

像这样声明矩阵，然后从那里开始：

M = np.zeros((n, n), dtype=float)
于 2019-12-14T22:44:21.973 回答