Python

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 不熟悉，而且我总体上缺乏编码技能，但我根本没有得到任何正确答案，我也不知道我哪里出错了。

我在这里添加了一个我正在使用的样本矩阵，并且这个 shoulddd 适用于，并且应该给出适当的 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

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

The input matrix A represents the superdiagonal half of the m×m Hermitian positive definite matrix to be factored.

所以需要替换输入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]]

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

The output matrix R represents the upper-triangular factor for which A = R' R

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

Python - 得到错误的 Cholesky 分解解？

Python - Getting the wrong solution for Cholesky decomposition?

math

matrix

linear-algebra