python – 稀疏矩阵的条件数

2023年1月4日 246次阅读

我正在尝试获取scipy稀疏矩阵的条件数.到目前为止我设法做到的方法是将矩阵转换为密集,然后获得它的特征值：

$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中正确完成这项工作？

最佳答案 如你所述,将稀疏矩阵转换为密集矩阵对于大问题通常不是一个好主意.因此,您不应该使用像numpy.linalg.cond()函数这样的函数来解决密集问题.

cond = || A || * || A ^ -1 ||

虽然它可能不是更有效的方法,但您可以在scipy.sparse的稀疏模块中使用函数normand inverse来评估条件数(反转矩阵是计算上昂贵的过程)：

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
点赞