python - 稀疏矩阵的条件数

Question

我正在尝试获取 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 中是如何正确完成的？

Answer 1

正如 Augusto Sisa 提到的，条件数定义为

条件(A):=||A|| ||A^{-1}||

但这很重要，即最大特征值的大小与最小(按大小)特征值的大小之比。因此，使用 scipy.sparse.linalg.eigs() Scipy reference manual 获取这些值更有意义并找出你自己。

import scipy.sparse
import scipy.sparse.linalg

ew1, ev = lg.eigsh(A, which='LM')
ew2, ev = lg.eigsh(K, sigma=1e-8)   #<--- takes a long time

ew1 = abs(ew1)
ew2 = abs(ew2)

condA = ew1.max()/ew2.min()

指定 sigma 选项是因为默认选项不能很好地找到最小的特征值。该选项以 shift-inverse 模式计算接近 sigma 的特征值。这一步需要很长时间。您可以通过将 k= 和 ncv= 选项指定为小于默认值的值来以牺牲准确性为代价来加快速度。例如 k=3 和 ncv=8。但是你真的不知道你是否会得到这些特征值的良好近似值。您几乎总是会通过使用 k 参数计算更多特征值来提高准确性，但对于大多数用途而言，默认值应该足够准确。亲自测试以查看矩阵的差异。

一般人们使用稀疏矩阵是因为矩阵的维度使得存储所有条目变得令人望而却步。如果数据存储对于给定的维度来说是令人望而却步的，那么对于一般矩阵来说，更糟糕的是逆计算。逆的计算本身是令人望而却步的，但是稀疏矩阵通常没有稀疏逆，所以你又遇到了第一个问题。