A = BT·B 的稀疏矩阵分解

Question

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

B.T.dot(B) = A

scipy模块上是否有任何功能可以做到这一点？如果没有，在numpy可能吗？

Answer 1

此功能在

numpy和scipy中可用于密集矩阵，scikits.sparse可用于稀疏矩阵，具体取决于您的矩阵。你的矩阵是正数吗？它是对称的吗？如果是这样，那么您本质上正在寻找乔列斯基分解。

>>> 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版本，您可以指定 keywork 参数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

类型分解：分析半定矩阵的Cholesky分解，但据我所知，没有在scipy或scikits中实现。