python - A = Bᵀ·B 的稀疏矩阵分解

Question

假设我有一个稀疏矩阵 A，我想计算一个矩阵 B 使得

B.T.dot(B) = A

scipy 模块中是否有任何函数可以做到这一点？如果不是，是否可以在 numpy 中实现？

Answer 1

此功能在 numpy 和 scipy 中适用于稠密矩阵，在 scikits.sparse 中适用于稀疏矩阵，具体取决于您的矩阵。您的矩阵正定吗？它是对称的吗？如果是这样，那么您实际上是在寻找 Cholesky 分解。

>>> A = np.random.rand(100,100)  # Construct a dense matrix
>>> np.fill_diagonal(A,10) # Ensure the matrix is positive definite
>>> A = 0.5*(A+A.T) # Symmetrize the matrix
>>> B = np.linalg.cholesky(A) # Perform Cholesky decomposition 
>>> np.allclose(B.dot(B.T),A) # Verify
True

通常 B 是 numpy 中的下三角矩阵，但最好确保这一点，至于上三角矩阵，您需要更改从 B.dot(B.T) 到 B.T.dot(B) 的乘法顺序。对于 scipy 版本，您可以指定关键字参数 lower(默认情况下为 False)，以始终获得正确的顺序需要。

>>> B = sp.linalg.cholesky(A,lower=False)
>>> np.allclose(B.T.dot(B),A)
True

虽然这应该没什么大不了的。对于稀疏矩阵，您可以使用 todense() 或 toarray()(不推荐)将稀疏矩阵转换为密集矩阵，或者使用 scikits.sparse模块

>>> from scikits.sparse.cholmod import cholesky
>>> spA = csc_matrix(A)
>>> factor = cholesky(spA)
>>> spB = factor.L()
>>> np.allclose(spA.todense(),spB.dot(spB.T).todense()) # Just for verification
True

对于半定矩阵，论文中也提供了 Cholesky 类型的分解:Analysis of the Cholesky Decomposition of a Semi-definite Matrix ，但据我所知没有在 scipy 或 scikits 中实现。