如何实现A=L^T·L形式的Cholesky分解？附实现代码

云上老周

2026-7-5

正定Hermite矩阵的Cholesky分解（A = L^H·L形式）

问题说明

numpy自带的np.linalg.cholesky函数返回的下三角矩阵L满足A = L·L^H**（LH表示L的共轭转置），但如果需要得到下三角矩阵L使得A = L^H·L**，就需要调整分解逻辑，自行实现对应的算法。

算法原理

针对正定Hermite矩阵的特性，我们可以通过逆序遍历矩阵行来计算下三角矩阵L的元素：

对角元素计算：从最后一行（i = n-1）开始往前遍历，L[i,i]的取值为：
$$L_{i,i} = \sqrt{A_{i,i} - \sum_{k=i+1}^{n-1} |L_{k,i}|^2}$$
非对角元素计算：对于i > j的位置，L[i,j]的取值为：
$$L_{i,j} = \frac{A_{i,j} - \sum_{k=i+1}^{n-1} \overline{L_{k,i}} \cdot L_{k,j}}{L_{i,i}}$$
其中$\overline{L_{k,i}}$表示$L_{k,i}$的共轭复数。

Python实现代码

import numpy as np

def cholesky_hermitian_adjoint(A):
    N = len(A)
    # 初始化复数类型的下三角矩阵
    L = np.zeros_like(A, dtype='complex_')
    
    # 逆序遍历行，从最后一行到第一行
    for i in reversed(range(N)):
        # 计算对角元素 L[i,i]
        sum_diag = 0.0
        for k in range(i+1, N):
            sum_diag += np.conjugate(L[k, i]) * L[k, i]
        L[i, i] = np.sqrt(A[i, i] - sum_diag)
        
        # 计算i行中j < i的非对角元素
        for j in reversed(range(i)):
            sum_offdiag = 0.0
            for k in range(i+1, N):
                sum_offdiag += np.conjugate(L[k, i]) * L[k, j]
            L[i, j] = (A[i, j] - sum_offdiag) / L[i, i]
    
    return L

# 示例使用
if __name__ == "__main__":
    # 3阶正定Hermite矩阵（实矩阵情况下等价于正定对称矩阵）
    A = np.array([[4, 12, -16],
                  [12, 37, -43],
                  [-16, -43, 98]], dtype=np.float64)
    
    L = cholesky_hermitian_adjoint(A)
    print("下三角矩阵L：")
    print(L)
    print("\n验证 L^H · L 是否等于原矩阵A：")
    print(np.matmul(np.conj(L.T), L))