Trouver les plus petits vecteurs propres de grande matrice clairsemée, plus de 100 fois plus lent dans SciPy que dans Octave

J'essaie de calculer quelques (5-500) vecteurs propres correspondant aux plus petites valeurs propres de grandes matrices creuses symétriques carrées (jusqu'à 30000x30000) avec moins de 0,1% des valeurs étant non nulles.

J'utilise actuellement scipy.sparse.linalg.eigsh en mode inversion de décalage (sigma = 0.0), ce que j'ai compris à travers divers articles sur le sujet est la solution préférée. Cependant, il faut jusqu'à 1h pour résoudre le problème dans la plupart des cas. D'un autre côté, la fonction est très rapide, si je demande les valeurs propres les plus élevées (sous-secondes sur mon système), ce qui était attendu de la documentation.

Comme je suis plus familier avec Matlab du travail, j'ai essayé de résoudre le problème dans Octave, ce qui m'a donné le même résultat en utilisant eigs (sigma = 0) en quelques secondes (moins de 10s). Étant donné que je veux faire un balayage des paramètres de l'algorithme, y compris le calcul du vecteur propre, ce type de gain de temps serait également formidable en python.

J'ai d'abord changé les paramètres (en particulier la tolérance), mais cela n'a pas beaucoup changé dans les délais.

J'utilise Anaconda sur Windows, mais j'ai essayé de changer le LAPACK / BLAS utilisé par scipy (ce qui était une énorme douleur) de mkl (Anaconda par défaut) à OpenBlas (utilisé par Octave selon la documentation), mais je n'ai pas pu voir de changement dans performance.

Je n'ai pas pu déterminer s'il y avait quelque chose à changer dans l'ARPACK utilisé (et comment)?

J'ai téléchargé un testcase pour le code ci-dessous dans le dossier déroulant suivant: https://www.dropbox.com/sh/l6aa6izufzyzqr3/AABqij95hZOvRpnnjRaETQmka?dl=0

En 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)

En octave:

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

Toute aide est appréciée!

Quelques options supplémentaires que j'ai essayées sur la base des commentaires et suggestions:

Octave: eigs(M,6,0)et eigs(M,6,'sm')donnez-moi le même résultat:

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

tout eigs(M,6,'sa',struct('tol',2))converge vers

[1.0423 2.7604 6.1548 11.1310 18.0207 25.3933] 

beaucoup plus rapide, mais uniquement si les valeurs de tolérance sont supérieures à 2, sinon elles ne convergent pas du tout et les valeurs sont fortement différentes.

Python: eigsh(M,k=6,which='SA')et les eigsh(M,k=6,which='SM')deux ne convergent pas (erreur ARPACK sur aucune convergence atteinte). eigsh(M,k=6,sigma=0.0)Donne seulement quelques valeurs propres (après presque une heure), qui sont différentes de l'octave pour les plus petites (même 1 petite valeur supplémentaire est trouvée):

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

Si la tolérance est suffisamment élevée, j'obtiens également des résultats eigsh(M,k=6,which='SA',tol='1')qui se rapprochent des autres valeurs obtenues

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

encore une fois avec un nombre différent de petites valeurs propres. Le temps de calcul est encore proche de 30min. Alors que les différentes valeurs très petites peuvent être compréhensibles, car elles peuvent représenter des multiples de 0, la multiplicité différente me déconcerte.

De plus, il semble y avoir des différences fondamentales dans SciPy et Octave, que je ne peux pas encore comprendre.

Une conjecture et quelques commentaires, puisque je n'ai pas Matlab / Octave:

Pour trouver de petites valeurs propres de matrices symétriques avec des valeurs propres> = 0, comme la vôtre, ce qui suit est waaay plus rapide que shift-invert:

# flip eigenvalues e.g.
# A:     0 0 0 ... 200 463
# Aflip: 0 163 ... 463 463 463
maxeval = eigsh( A, k=1 )[0]  # biggest, fast
Aflip = maxeval * sparse.eye(n) - A
bigevals, evecs = eigsh( Aflip, which="LM", sigma=None ... )  # biggest, near 463
evals = maxeval - bigevals  # flip back, near 463 -> near 0
# evecs are the same

eigsh( Aflip )pour les grandes paires propres, c'est plus rapide que shift-invert pour les petits, parce que A * xc'est plus rapide que solve()ce shift-invert doit faire. Matlab / Octave pourrait éventuellement le faire Aflipautomatiquement, après un test rapide de positif-défini avec Cholesky.
Pouvez-vous exécuter eigsh( Aflip )dans Matlab / Octave?

Autres facteurs pouvant affecter la précision / vitesse:

La valeur par défaut d'Arpack pour le vecteur de départ v0est un vecteur aléatoire. J'utilise v0 = np.ones(n), ce qui peut être terrible pour certains Amais est reproductible :)

Cette Amatrice est presque exactement sigulaire, A * ones~ 0.

Multicore: scipy-arpack avec openblas / Lapack utilise environ 3,9 des 4 cœurs de mon iMac; Matlab / Octave utilise-t-il tous les cœurs?

k 10  tol 1e-05:    8 sec  eigvals [0 8.5e-05 0.00043 0.0014 0.0026 0.0047 0.0071 0.0097 0.013 0.018] 
k 10  tol 1e-06:   44 sec  eigvals [0 3.4e-06 2.8e-05 8.1e-05 0.00015 0.00025 0.00044 0.00058 0.00079 0.0011] 
k 10  tol 1e-07:  348 sec  eigvals [0 3e-10 7.5e-07 7.6e-06 1.2e-05 1.9e-05 2.1e-05 4.2e-05 5.7e-05 6.4e-05] 

k 20  tol 1e-05:   18 sec  eigvals [0 5.1e-06 4.5e-05 0.00014 0.00023 0.00042 0.00056 0.00079 0.0011 0.0015 0.0017 0.0021 0.0026 0.003 0.0037 0.0042 0.0047 0.0054 0.006
k 20  tol 1e-06:   73 sec  eigvals [0 5.5e-07 7.4e-06 2e-05 3.5e-05 5.1e-05 6.8e-05 0.00011 0.00014 0.00016 0.0002 0.00025 0.00027 0.0004 0.00045 0.00051 0.00057 0.00066
k 20  tol 1e-07:  267 sec  eigvals [-4.8e-11 0 7.5e-07 7.6e-06 1e-05 1.9e-05 2e-05 2.2e-05 4.2e-05 5.1e-05 5.8e-05 6.4e-05 6.9e-05 8.3e-05 0.00011 0.00012 0.00013 0.00015

k 50  tol 1e-05:   82 sec  eigvals [-4e-13 9.7e-07 1e-05 2.8e-05 5.9e-05 0.00011 0.00015 0.00019 0.00026 0.00039 ... 0.0079 0.0083 0.0087 0.0092 0.0096 0.01 0.011 0.011 0.012
k 50  tol 1e-06:  432 sec  eigvals [-1.4e-11 -4e-13 7.5e-07 7.6e-06 1e-05 1.9e-05 2e-05 2.2e-05 4.2e-05 5.1e-05 ... 0.00081 0.00087 0.00089 0.00096 0.001 0.001 0.0011 0.0011
k 50  tol 1e-07: 3711 sec  eigvals [-5.2e-10 -4e-13 7.5e-07 7.6e-06 1e-05 1.9e-05 2e-05 2.2e-05 4.2e-05 5.1e-05 ... 0.00058 0.0006 0.00063 0.00066 0.00069 0.00071 0.00075

versions: numpy 1.18.1  scipy 1.4.1  umfpack 0.3.2  python 3.7.6  mac 10.10.5 

Matlab / Octave sont-ils à peu près les mêmes? Sinon, tous les paris sont désactivés - vérifiez d'abord l'exactitude, puis la vitesse.

Pourquoi les valeurs propres vacillent-elles autant? Un petit <0 pour une matrice supposée non négative-définie est un signe d' erreur d'arrondi , mais l'astuce habituelle d'un minuscule décalage A += n * eps * sparse.eye(n), n'aide pas.

J'espère que cela t'aides.

Je sais que c'est vieux maintenant, mais j'ai eu le même problème. Avez-vous révisé ici ( https://docs.scipy.org/doc/scipy/reference/tutorial/arpack.html )?

Il semble que lorsque vous définissez sigma sur un nombre faible (0), vous devez définir which = 'LM', même si vous souhaitez des valeurs faibles. En effet, la définition de sigma transforme les valeurs que vous souhaitez (faibles dans ce cas) en apparence élevées et vous pouvez donc toujours profiter des méthodes `` LM '', qui sont beaucoup plus rapides pour obtenir ce que vous voulez (les valeurs propres faibles ).

Je veux d'abord dire que je ne sais pas pourquoi les résultats que vous et @Bill avez rapportés sont tels qu'ils sont. Je me demande simplement si eigs(M,6,0)dans Octave correspond à k=6 & sigma=0, ou peut-être que c'est autre chose?

Avec scipy, si je ne fournis pas de sigma, je peux obtenir un résultat en un temps décent de cette façon.

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)

Je ne suis absolument pas sûr que cela ait du sens cependant.

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]

Le seul moyen que j'ai trouvé d'utiliser sigma et d'obtenir un résultat dans un temps décent est de fournir M comme opérateur linéaire. Je ne suis pas trop familier avec cette chose, mais d'après ce que j'ai compris, mon implémentation représente une matrice d'identité, qui est ce que M devrait être s'il n'est pas spécifié dans l'appel. La raison en est qu'au lieu d'effectuer une résolution directe (décomposition LU), scipy utilisera un solveur itératif, qui est potentiellement mieux adapté. À titre de comparaison, si vous fournissez M = np.identity(a.shape[0])ce qui devrait être exactement le même, eigsh prend une éternité pour produire un résultat. Notez que cette approche ne fonctionne pas si elle sigma=0est fournie. Mais je ne sais pas si sigma=0est-ce vraiment utile?

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)))

Encore une fois, aucune idée si elle est correcte mais définitivement différente d'avant. Ce serait formidable d'avoir la contribution de quelqu'un de 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]