找到大型稀疏矩阵的最小特征向量,在SciPy中比在八度中慢100倍以上

问题描述 投票:7回答:1

我正在尝试计算与大型对称方形稀疏矩阵(最大30000x30000)的最小特征值相对应的很少的(5-500)特征向量,其中小于0.1%的特征值为非零。

我目前正在以移位-反转模式(sigma = 0.0)使用scipy.sparse.linalg.eigsh,我通过有关该主题的各种帖子发现这是首选的解决方案。但是,在大多数情况下,最多需要1小时才能解决该问题。另一方面,如果我要求从文档中获得最大的特征值(系统中的秒数),则该功能将非常快。

由于我对工作更熟悉Matlab,所以我尝试在Octave中解决问题,这在短短几秒钟(不到10秒)内使用eigs(sigma = 0)给出了相同的结果。由于我想对包括特征向量计算在内的算法进行参数扫描,因此在python中也要具有这种时间增益。

我首先更改了参数(尤其是公差),但是在时间尺度上并没有太大变化。

我正在Windows上使用Anaconda,但尝试将scipy(这是一个巨大的痛苦)使用的LAPACK / BLAS从mkl(默认为Anaconda)切换为OpenBlas(根据文档由Octave使用),但看不到性能变化。

我不知道使用过的ARPACK是否有某些更改(以及如何更改?)>

我将以下代码的测试用例上传到了以下保管箱文件夹:https://www.dropbox.com/sh/l6aa6izufzyzqr3/AABqij95hZOvRpnnjRaETQmka?dl=0

在Python中

import numpy as np
from scipy.sparse import csr_matrix, csc_matrix, linalg, load_npz   
M = load_npz('M.npz')
evals, evecs = linalg.eigsh(M,k=6,sigma=0.0)

在八度音阶中:

M=dlmread('M.txt');
M=spconvert(M);
[evecs,evals] = eigs(M,6,0);

需要任何帮助!

我根据评论和建议尝试了一些其他选项:

八度:eigs(M,6,0)eigs(M,6,'sm')给我相同的结果:

[1.8725e-05 1.0189e-05 7.5622e-06 7.5420e-07 -1.2239e-18 -2.5674e-16]

[eigs(M,6,'sa',struct('tol',2))收敛到]时>

[1.0423 2.7604 6.1548 11.1310 18.0207 25.3933] 

要快得多,但是仅当公差值大于2时,否则它根本不会收敛,并且值有很大差异。

Python:eigsh(M,k=6,which='SA')eigsh(M,k=6,which='SM')都不收敛(在未收敛时出现ARPACK错误)。只有eigsh(M,k=6,sigma=0.0)给出一些特征值(近一个小时后),这些特征值与最小特征值的倍频程不同(甚至找到了另外一个小值):

[3.82923317e-17 3.32269886e-16 2.78039665e-10 7.54202273e-07 7.56251500e-06 1.01893934e-05]

如果公差足够高,我也会从eigsh(M,k=6,which='SA',tol='1')中获得结果,接近其他获得的值

[4.28732218e-14 7.54194948e-07 7.56220703e-06 1.01889544e-05, 1.87247350e-05 2.02652719e-05]

再次具有不同数量的小特征值。计算时间仍然接近30分钟。尽管可以理解不同的非常小的值,但由于它们可能表示0的倍数,因此不同的多重性使我感到困惑。

此外,SciPy和Octave之间似乎存在一些根本差异,但我无法弄清楚。

我正在尝试计算少量(5-500)特征向量,这些特征向量对应于大型对称方形稀疏矩阵(最大30000x30000)的最小特征值,其中小于0.1%的特征值非零。...]] >

我首先要说的是,我不知道您和@Bill报告的结果为何如此。我只是想知道八度音阶中的eigs(M,6,0)是否对应于k=6 & sigma=0,还是其他?

使用scipy,如果我不提供sigma,则可以通过这种方式在适当的时间获得结果。

import numpy as np
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import eigsh
from time import perf_counter
M = np.load('M.npz')
a = csr_matrix((M['data'], M['indices'], M['indptr']), shape=M['shape'])
t = perf_counter()
b, c = eigsh(a, k=50, which='SA', tol=1e1)
print(perf_counter() - t)
print(b)

我完全不确定这是否有意义。

0.4332823531003669
[4.99011753e-03 3.32467891e-02 8.81752215e-02 1.70463893e-01
 2.80811313e-01 4.14752072e-01 5.71103821e-01 7.53593653e-01
 9.79938915e-01 1.14003837e+00 1.40442848e+00 1.66899183e+00
 1.96461415e+00 2.29252666e+00 2.63050114e+00 2.98443218e+00
 3.38439528e+00 3.81181747e+00 4.26309942e+00 4.69832271e+00
 5.22864462e+00 5.74498014e+00 6.22743988e+00 6.83904055e+00
 7.42379697e+00 7.97206446e+00 8.62281827e+00 9.26615266e+00
 9.85483434e+00 1.05915030e+01 1.11986296e+01 1.18934953e+01
 1.26811461e+01 1.33727614e+01 1.41794599e+01 1.47585155e+01
 1.55702295e+01 1.63066947e+01 1.71564622e+01 1.78260727e+01
 1.85693454e+01 1.95125277e+01 2.01847294e+01 2.09302671e+01
 2.18860389e+01 2.25424795e+01 2.32907153e+01 2.37425085e+01
 2.50784800e+01 2.55119112e+01]

我发现使用sigma并在适当的时间内获得结果的唯一方法是将M提供为LinearOperator。我对此不太熟悉,但是据我了解,我的实现表示一个单位矩阵,如果未在调用中指定M,则应为M。这样做的原因是scipy不会执行直接求解(LU分解),而是使用迭代求解器,它可能更适合。作为比较,如果您提供的M = np.identity(a.shape[0])应该完全相同,则eigsh会花费很多时间才能得出结果。请注意,如果提供sigma=0,则此方法无效。但是我不确定sigma=0是否真的有用吗?

import numpy as np
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import eigs, eigsh, LinearOperator
from time import perf_counter


def mv(v):
    return v


M = np.load('M.npz')
a = csr_matrix((M['data'], M['indices'], M['indptr']), shape=M['shape'])
t = perf_counter()
b, c = eigsh(a, M=LinearOperator(shape=a.shape, matvec=mv, dtype=np.float64),
             sigma=5, k=50, which='SA', tol=1e1, mode='cayley')
print(perf_counter() - t)
print(np.sort(-5 * (1 + b) / (1 - b)))

再次,不知道它是否正确,但肯定与以前有所不同。希望有人得到scipy的帮助。

1.4079377939924598
[3.34420263 3.47938816 3.53019328 3.57981026 3.60457277 3.63996294
 3.66791416 3.68391585 3.69223712 3.7082205  3.7496456  3.76170023
 3.76923989 3.80811939 3.81337342 3.82848729 3.84137264 3.85648208
 3.88110869 3.91286153 3.9271108  3.94444577 3.97580798 3.98868207
 4.01677424 4.04341426 4.05915855 4.08910692 4.12238969 4.15283192
 4.16871081 4.1990492  4.21792125 4.24509036 4.26892806 4.29603036
 4.32282475 4.35839271 4.37934257 4.40343219 4.42782208 4.4477206
 4.47635849 4.51594603 4.54294049 4.56689989 4.58804775 4.59919363
 4.63700551 4.66638214]
python scipy octave sparse-matrix eigenvector
1个回答
0
投票

我首先要说的是,我不知道您和@Bill报告的结果为何如此。我只是想知道八度音阶中的eigs(M,6,0)是否对应于k=6 & sigma=0,还是其他?

使用scipy,如果我不提供sigma,则可以通过这种方式在适当的时间获得结果。

© www.soinside.com 2019 - 2024. All rights reserved.