稀疏矩阵的条件数

Question

我正在尝试获取稀疏矩阵的条件编号。到目前为止，我设法做到的方法是将矩阵转换为稠密的，然后获得其特征值：

$ python
Python 3.5.2 (v3.5.2:4def2a2901a5, Jun 26 2016, 10:47:25) 
[GCC 4.2.1 (Apple Inc. build 5666) (dot 3)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> from numpy import array
>>> import numpy as np
>>> import scipy.sparse as sparse
>>> I = array([0,3,1,0])
>>> J = array([0,3,1,2])
>>> V = array([4,5,7,9])
>>> A = sparse.coo_matrix((V,(I,J)),shape=(4,4))
>>> A = A.todense()
>>> eig = np.linalg.eig(A)
>>> eig = eig[0].real, np.array(eig[1].real)
>>> def split(array, cond):
...     return (array[cond], array[~cond])
... 
>>> eigv, zero = split(eig[0], eig[0]>1e-10)
>>> cond = max(eigv) / min(eigv)
>>> cond
1.75

如您所料，这对于大型矩阵将变得不可行。我想知道如何在Python中正确完成此操作？

Answer 1

正如您所提到的，对于大型问题，将稀疏矩阵转换为密集矩阵通常不是一个好主意。因此，您不应该使用类似numpy.linalg.cond()对密集问题有用的函数之类的函数。

条件= || 一个 || * || A^-1 ||

尽管这可能不是更有效的方法，但您可以在稀疏模块中使用函数normand来评估条件数（反转矩阵是一个计算成本高昂的过程）：inversescipy.sparse

norm_A = scipy.sparse.linalg.norm(A)
norm_invA = scipy.sparse.linalg.norm(scipy.sparse.linalg.inv(A))
cond = norm_A*norm_invA

例子：

import numpy as np
import scipy.sparse as sparse

n = 7
diagonals = np.array([[1, -4, 6, -4, 1]]).repeat(n, axis=0).transpose()
offsets = np.array([-2, -1, 0, 1, 2])
S = sparse.dia_matrix((diagonals, offsets), shape=(n, n))
norm_S = sparse.linalg.norm(S, np.inf)
norm_invS = sparse.linalg.norm(sparse.linalg.inv(S), np.inf)
cond = norm_S*norm_invS