使用一些 Python 库对正半定矩阵进行 Cholesky 分解

Question

我正在寻找一个内置函数，可以在 Python 中实现正半定矩阵的 Cholesky 分解。

存在 NumPy (

numpy.linalg.cholesky

) 和 SciPy (

scipy.linalg.cholesky

) 的实现，但它们仅适用于正定矩阵。

我正在寻找正半定矩阵的扩展，但没有找到。

Answer 1

您可以使用

scipy.linalg.lapack.dpstrf

，它是 LAPACK

dpstrf

的漂亮骨架包装。

import numpy as np
import scipy.linalg

# Generate a positive semidefinite matrix
n, m = 10, 7
A = np.random.rand(n, m)
A = A @ A.T

L, piv, rank, info = scipy.linalg.lapack.dpstrf(A, lower=1)

# For clarity, generate explicit permutation matrix based on definition in [2]
# (For efficiency, we would permute the matrix directly)
P = np.zeros((n, n))
P[piv-1, np.arange(n)] = 1

# We've set `lower=1`; by default it would be upper-triangular
L = np.tril(L)

# Information about complete-pivot Cholesky suggests that
# we should end up with all-zero columns, but that's not
# what we see in LAPACK's output. Apparently this is needed
# if you want to check the decomposition below.
L[:, rank:] = 0

np.testing.assert_equal(rank, m)
np.testing.assert_allclose(P.T @ A @ P, L @ L.T)