Python - 为 Cholesky 分解找到错误的解决方案？

Question

我正在尝试将一些伪代码从 matlab 转换为 python 脚本，但在获得正确答案时遇到了一些麻烦？谁能帮我确定我在哪里搞砸了翻译？

我在 Trefethan & Bau 书中给出的用于 Cholesky 分解的伪代码是

如果我正确理解了给出的描述，这是为上三角矩阵完成的，但我认为这也适用于一般矩阵，不是吗？

反正我用python写了如下代码：

def is_SPD(A):
    if np.all(A == A.T):
        if np.all(la.eigvals(A) > 0):
            return True
    return False

def cholesky_decomp(A):
    #first check if matrix is symmetric and positive definite
    if is_SPD(A) == True:
        R = np.copy(A)
        for k in range (0, len(A)):
            for j in range (k+1, len(A)):
                R[j,j:] = R[j,j:] - (R[k,j:]*(R[k,j]/R[k,k]))
            R[k,k:] = R[k,k:]/sqrt(R[k,k])
        return R
    else:
        return print('Cholesky decomposition not applicable')

我在做一个4x4的矩阵，我用np.linalg方法做了分解，我的答案完全不一样。

我认为这可能吗？是因为我不熟悉 MATLAB，而且我通常缺乏编码技能，但我根本没有得到任何正确的答案，而且我不知道我哪里出错了。

我在这里添加了一个我正在使用的示例矩阵，并且这个 应该 适用于，并且应该给出适当的 Cholesky 分解，但我得到了一个完全错误的答案.

有人可以用它来帮助我找出哪里出错了吗？

A = np.array([[16, -12, -12, -16], [-12, 25, 1, -4], [-12, 1, 17, 14], [-16, -4, 14, 57]])

我的代码给了我：

[[  4  -3  -3  -4]
 [-12   4  -2  -4]
 [-12   1   2  -3]
 [-16  -4  14   4]]

虽然 numpy Cholesky 函数给了我：

[[ 4.  0.  0.  0.]
 [-3.  4.  0.  0.]
 [-3. -2.  2.  0.]
 [-4. -4. -3.  4.]]

Answer 1

您的代码正确地实现了规定的算法，但请注意文本显示（添加强调）：

输入矩阵A表示m×m厄米正定矩阵的上对角线一半为因素。

所以需要替换输入的A，

[[ 16 -12 -12 -16]
 [-12  25   1  -4]
 [-12   1  17  14]
 [-16  -4  14  57]]

np.triu(A):

[[ 16 -12 -12 -16]
 [  0  25   1  -4]
 [  0   0  17  14]
 [  0   0   0  57]]

此外（符号略有改变，其中' 表示厄米转置），

输出矩阵R代表上三角因子，其中A = R' R

所以你得到一个上三角 R，而 Numpy 的 cholesky 函数给出一个下三角结果。这些结果中的每一个都只是另一个的 Hermitian 转置版本（请参阅here）。