使用scipy求逆大稀疏矩阵

Question

我需要对一个很大的稀疏矩阵进行求逆。这个矩阵的求逆是必须的，唯一的捷径就是只关注主对角线上的元素，忽略其他的元素（虽然我不太想这样做，但如果能解决问题也可以接受）。

我需要求逆的矩阵通常很大（比如40000 * 40000），而且只有少数几个非零的对角线元素。我现在的做法是把所有东西都做成稀疏矩阵，然后

posterior_covar = np.linalg.inv ( hessian.todense() )

这样做显然需要很长时间，还占用很多内存。

有没有什么建议，还是说这只是需要耐心，或者把问题缩小一点呢？

Answer 1

我觉得稀疏模块没有明确的逆操作方法，但它确实有稀疏求解器。像这个简单的例子就能正常工作：

>>> a = np.random.rand(3, 3)
>>> a
array([[ 0.31837307,  0.11282832,  0.70878689],
       [ 0.32481098,  0.94713997,  0.5034967 ],
       [ 0.391264  ,  0.58149983,  0.34353628]])
>>> np.linalg.inv(a)
array([[-0.29964242, -3.43275347,  5.64936743],
       [-0.78524966,  1.54400931, -0.64281108],
       [ 1.67045482,  1.29614174, -2.43525829]])

>>> a_sps = scipy.sparse.csc_matrix(a)
>>> lu_obj = scipy.sparse.linalg.splu(a_sps)
>>> lu_obj.solve(np.eye(3))
array([[-0.29964242, -0.78524966,  1.67045482],
       [-3.43275347,  1.54400931,  1.29614174],
       [ 5.64936743, -0.64281108, -2.43525829]])

注意，结果是转置的！

如果你希望得到的逆矩阵也是稀疏的，而最后的求解结果太大，无法放进内存里，你还可以选择逐行（或逐列）生成它，提取出非零值，然后根据这些值来构建稀疏的逆矩阵：

>>> for k in xrange(3) :
...     b = np.zeros((3,))
...     b[k] = 1
...     print lu_obj.solve(b)
... 
[-0.29964242 -0.78524966  1.67045482]
[-3.43275347  1.54400931  1.29614174]
[ 5.64936743 -0.64281108 -2.43525829]