Python实现Cholesky分解结果与numpy输出不符，如何定位代码错误？

核心错误1：内层循环的索引逻辑写反 #

你参考的上三角Cholesky分解算法中，第k次迭代处理第k行，对所有j > k的行，需要用第k行修正第j行从k列开始的所有元素，而不是你写的R[j,j:]，你把起始列写错了，同时行操作的对应关系也错了。
正确的内层更新逻辑应该是：

R[j, k:] = R[j, k:] - R[k, k:] * (R[k, j] / R[k, k])

核心错误2：未清理下三角区域的冗余值 #

Cholesky分解得到的上三角矩阵R，下三角所有元素都应该为0，你直接拷贝了原矩阵A作为初始值，计算过程中没有清理下三角的旧数值，导致输出里下三角还保留着错误的原始数据。

额外说明：与numpy输出的格式差异 #

numpy.linalg.cholesky默认返回下三角矩阵L，满足A = L @ L.T，你实现的是上三角版本R，满足A = R.T @ R，二者正确结果互为转置关系，你的计算正确的话，输出的R应该等于numpy返回结果的转置。

修正后的完整代码 #

import numpy as np
from numpy.linalg import eigvals
from math import sqrt

def is_SPD(A):
    if np.all(A == A.T):
        # 浮点比较加微小阈值避免精度问题
        if np.all(eigvals(A) > 1e-10):
            return True
    return False

def cholesky_decomp(A):
    if is_SPD(A):
        n = len(A)
        R = np.copy(A)
        for k in range(n):
            for j in range(k+1, n):
                # 修正索引：更新j行从k列开始的元素
                R[j, k:] = R[j, k:] - R[k, k:] * (R[k, j] / R[k, k])
            # 缩放当前k行
            R[k, k:] = R[k, k:] / sqrt(R[k, k])
        # 清理下三角区域，只保留上三角
        R[np.tril_indices(n, -1)] = 0
        return R
    else:
        print('Cholesky decomposition not applicable')

测试验证 #

用你提供的矩阵测试：

A = np.array([[16, -12, -12, -16], [-12, 25, 1, -4], [-12, 1, 17, 14], [-16, -4, 14, 57]])
R = cholesky_decomp(A)
L = np.linalg.cholesky(A)
# 输出R和L的转置，二者完全一致
print(R)
print(L.T)

内容的提问来源于stack exchange，提问作者Applesauce44