具有互换行和列的稀疏 Cholesky 分解答案

【问题标题】：sparse cholesky decomposition with interchanged rows and columns具有互换行和列的稀疏 Cholesky 分解
【发布时间】：2016-07-10 20:24:00
【问题描述】：

我正在使用 python 的 scikits.sparse.cholmod 来获得对称矩阵的 Cholesky 分解。

我将 cholesky() 的结果与 matlab 的 chol() 进行了比较。结果与某些行和列互换时存在差异。我正在尝试迭代分解以获得特征值，这种差异似乎是有问题的。

这是我的代码：

import numpy as np
from scipy.sparse import csr_matrix
from scipy.sparse import csc_matrix
from scikits.sparse.cholmod import cholesky

A = csr_matrix([[1,2,0,0], [0,0,3, 0], [4,0,5, 0], [0, 0, 1, 2]])
B = (A*A.T)
print "B: "
print B.todense()

for i in range(10):
    factor = cholesky(B.tocsc())
    l = factor.L()  #l is lower triangular
    B = (l.T*l)
    print l.todense()

第一次迭代的下三角矩阵是：

[[ 2.23606798  0.         0.          0.        ]
[ 0.          3.          0.          0.        ]
[ 0.          1.          2.          0.        ]
[ 1.78885438  5.          0.          3.57770876]]

而matlab的下三角矩阵为：

[2.2361        0         0         0
     0    3.0000         0         0
1.7889    5.0000    3.5777         0
     0    1.0000         0    2.0000]

matlab 结果是合理的，因为它导致正确的特征值。我在 python 中选择稀疏矩阵的类型有问题吗？

【问题讨论】：

标签： python sparse-matrix matrix-factorization

【解决方案1】：

cholesky 算法正在使用填充减少算法。因此，它设置了一个置换矩阵P。这样LL'=PBP'。

您可以参考factor documentation了解更多信息。

如果你打印P，你会得到：

>>> factor.P()
array([0, 1, 3, 2], dtype=int32)

这正是两个矩阵之间的差异。最后两行和最后两列的排列。

【讨论】：