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

12

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.

Spacekiller23
la source
2
1 - Je suppose que vous vouliez mettre des crochets autour de [evals, evecs] dans le code d'octave? 2 - pouvez-vous inclure un petit exemple pour M? ou peut-être un script générateur pour un si c'est possible?
Nick J
1 - Oui exactement, j'ai édité mon post. 2 - J'ai essayé les performances de certaines sous-matrices de mes données et il semble qu'Octave soit toujours plus rapide, mais pour les petites matrices inférieures à 5000x5000, c'est juste un facteur de 2 à 5 fois, au-dessus, cela devient vraiment moche. Et depuis ses "vraies données" je ne peux pas donner de script générateur. Existe-t-il un moyen standard de télécharger un exemple d'une manière ou d'une autre? En raison de la rareté, un fichier npz est raisonnablement petit.
Spacekiller23
Je suppose que vous pouvez partager un lien vers n'importe quelle installation de stockage en nuage.
Patol75
THX. J'ai inclus un lien dropbox dans le message d'origine et mis à jour le code en un exemple de travail.
Spacekiller23
1
Pour renforcer votre argument, j'ai testé sur Matlab R2019b et obtenu 84 secondes contre 36,5 minutes en Python 3.7, Scipy 1.2.1 (26 fois plus rapide).
Bill

Réponses:

1

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?


Voici les valeurs propres de scipy-Arpack pour plusieurs ket tol, récupérées à partir des fichiers journaux sous gist.github :

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.


D'où cela Avient-il, quel problème? Pouvez-vous générer des Aimages similaires , plus petites ou plus rares?

J'espère que cela t'aides.

denis
la source
Merci pour votre contribution et désolé pour la réponse (très) tardive. Le projet pour lequel j'ai utilisé ceci est déjà terminé, mais je suis toujours curieux, alors j'ai vérifié. Malheureusement, les valeurs propres dans Ocatve s'avèrent différentes, pour k = 10, je trouve [-2.5673e-16 -1.2239e-18 7.5420e-07 7.5622e-06 1.0189e-05 1.8725e-05 2.0265e-05 2.1568e- 05 4.2458e-05 5.1030e-05] qui est également indépendante de la valeur de tolérance dans la plage 1e-5 à 1e-7. Il y a donc une autre différence ici. Ne pensez-vous pas qu'il est étrange que scipy (y compris votre suggestion) produise différentes petites valeurs en fonction du nombre de valeurs interrogées?
Spacekiller23
@ Spacekiller23, c'était un bug, maintenant corrigé dans scipy 1.4.1 (voir scipy / issues / 11198 ); pourriez-vous vérifier votre version? Est également toldésordonné pour les petites valeurs propres - posez une nouvelle question à ce sujet si vous le souhaitez, faites le moi savoir.
denis
1

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

Anthony Gatti
la source
Cela a-t-il réellement changé les performances pour vous? Ce serait une surprise pour moi. Je connaissais le lien que vous publiez et j'ai également implicitement spécifié lequel = 'LM' dans mon exemple. Parce que la valeur par défaut pour un unset qui est 'LM'. J'ai quand même vérifié, et les performances sont inchangées pour mon exemple.
Spacekiller23
En effet, semble avoir une différence similaire à celle de Python à l'octave pour vous. J'avais également une grande matrice que j'essayais de décomposer et j'ai fini par utiliser eigsh (matrice, k = 7, qui = 'LM', sigma = 1e-10). À l'origine, je spécifiais à tort quel = 'SM' pensais que je devais faire cela pour obtenir les plus petites valeurs propres et cela prenait une éternité pour me donner une réponse. Ensuite, j'ai trouvé cet article et j'ai réalisé que vous aviez juste besoin de le spécifier sur le «LM» le plus rapide, et de définir k pour qu'il soit ce que vous vouliez et cela accélérerait les choses. Votre matrice est-elle réellement hermitienne?
Anthony Gatti
0

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]
Patol75
la source
Merci pour votre contribution et vos commentaires. J'ai essayé certaines choses pour donner une réponse décente à vos points. 1. Ma tâche consiste à trouver les k plus petites valeurs / vecteurs propres. Par conséquent, l'approche utilisant sigma = 0 est même donnée dans les documents SciPy: docs.scipy.org/doc/scipy/reference/tutorial/arpack.html 2. J'ai essayé quelques autres options, que j'ai éditées dans la question d'origine. 3. Si je comprends bien les documentaires d'Octave et SciPy, eigs (M, 6,0) et k = 6, simga = 0 devrait être le même.
Spacekiller23
4. Puisque ma matrice est réelle et carrée, je pensais qu'il ne devrait pas y avoir de différence entre SA et SM en option, mais il y en a évidemment, du moins dans le calcul. Suis-je sur une mauvaise piste ici? Dans l'ensemble, cela signifie plus de questions et mais pas de vraies réponses ou solutions de ma part.
Spacekiller23