我正在使用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中选择稀疏矩阵的类型是不是做错了什么?
cholesky算法使用的是填充缩减算法。因此,它建立了一个置换矩阵
有关更多信息,您可以参考因子文档。
如果打印
>>> factor.P()
array([0, 1, 3, 2], dtype=int32)
这正是两个矩阵之间的差异。最后两行和两列的排列。