Scipy/Numpy：执行Cholesky分解时检查矩阵正定性

嘿，这个问题我之前把Matlab脚本转Python的时候也踩过坑！Matlab的chol函数确实很贴心，能同时返回分解结果和矩阵的正定状态，不会一遇到浮点近似误差导致的非正定就直接报错，但Scipy默认的Cholesky分解确实会抛出异常。不过别担心，我们可以用几种方式实现类似Matlab的容错行为：

方法1：基础容错版（返回正定标记） #

这个方法最贴近Matlab [L,p] = chol(A) 的核心逻辑——先尝试原矩阵的分解，失败后就添加对角扰动再试，同时返回矩阵是否原本就是正定的标记：

import numpy as np
from scipy.linalg import cholesky, LinAlgError

def chol_with_positive_check(mat, eps=1e-8):
    try:
        # 尝试对原矩阵做Cholesky分解
        L = cholesky(mat, lower=True)
        return L, True  # 返回分解矩阵 + 原矩阵正定的标记
    except LinAlgError:
        # 原矩阵非正定，尝试添加小对角扰动
        mat_perturbed = mat + np.eye(mat.shape[0]) * eps
        try:
            L = cholesky(mat_perturbed, lower=True)
            return L, False  # 返回扰动后的分解矩阵 + 原矩阵非正定的标记
        except LinAlgError:
            # 加扰动后仍无法分解，说明矩阵本身存在严重问题
            raise ValueError("矩阵在添加对角扰动后仍无法完成Cholesky分解")

使用示例：

# 构造一个接近正定的矩阵（模拟浮点误差导致的非正定）
A = np.array([[1.0, 0.9999999], [0.9999999, 1.0]])
L, is_positive_definite = chol_with_positive_check(A)
print(f"原矩阵是否正定？{is_positive_definite}")

方法2：复刻Matlab的p值返回（定位非正定主子式） #

如果你想完全复刻Matlab中[L,p] = chol(A)的行为——p=0表示矩阵正定，p>0表示第p阶主子式非正定，可以用下面的函数：

import numpy as np
from scipy.linalg import cholesky, LinAlgError

def chol_with_p_value(mat, eps=1e-8):
    n = mat.shape[0]
    # 逐个检查主子式
    for p in range(1, n+1):
        submat = mat[:p, :p]
        try:
            cholesky(submat, lower=True)
        except LinAlgError:
            # 当前主子式非正定，尝试添加扰动
            submat_perturbed = submat + np.eye(p) * eps
            try:
                cholesky(submat_perturbed, lower=True)
                # 对整个矩阵添加扰动后完成分解
                L = cholesky(mat + np.eye(n)*eps, lower=True)
                return L, p
            except LinAlgError:
                raise ValueError(f"第{p}阶主子式在添加扰动后仍无法完成Cholesky分解")
    # 所有主子式都正定，正常分解
    L = cholesky(mat, lower=True)
    return L, 0

这个函数返回的p值和Matlab完全对应，方便你定位问题出在矩阵的哪个部分。

注意事项 #

• 扰动的eps值可以根据你的矩阵规模和精度需求调整，一般1e-8到1e-6是比较常用的范围；
• 如果矩阵本身就是严重非正定（不是浮点误差导致的），那添加扰动也无法解决问题，这时候抛出异常是合理的，避免后续计算出现错误；
• 如果你只需要分解结果而不关心原矩阵是否正定，也可以直接使用scipy.linalg.cholesky并配合try-except块处理异常，省去返回标记的逻辑。

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