Quelle est la façon la plus rapide de calculer la plus grande valeur propre d'une matrice générale?

27

EDIT: Je teste si des valeurs propres ont une magnitude de un ou plus.

J'ai besoin de trouver la plus grande valeur propre absolue d'une grande matrice clairsemée et non symétrique.

J'ai utilisé la eigen()fonction de R , qui utilise l'algo QR d'EISPACK ou de LAPACK pour trouver toutes les valeurs propres, puis j'utilise abs()pour obtenir les valeurs absolues. Cependant, je dois le faire plus rapidement.

J'ai également essayé d'utiliser l'interface ARPACK dans le igraphpackage R. Cependant, cela a donné une erreur pour l'une de mes matrices.

L'implémentation finale doit être accessible depuis R.

Il y aura probablement plusieurs valeurs propres de même ampleur.

Avez-vous des suggestions?

EDIT: La précision doit seulement être 1e-11. Une matrice "typique" a jusqu'à présent été . J'ai pu y faire une factorisation QR. Cependant, il est également possible d'en avoir de plus grandes. Je commence actuellement à lire sur l'algorithme Arnoldi. Je comprends que cela est lié à Lanczsos. $386\times 386$

EDIT2: Si j'ai plusieurs matrices que je "teste" et que je sais qu'il existe une grande sous-matrice qui ne varie pas. Est-il possible de l'ignorer / de la rejeter?

linear-algebra algorithms eigensystem sparse-matrix Puissance
la source

Voir ma réponse ici: scicomp.stackexchange.com/a/1679/979 . Il s'agit d'un sujet de recherche actuel et les méthodes actuelles peuvent faire mieux que Lanczos. Le problème du calcul des valeurs singulières est équivalent au problème du calcul des valeurs propres.

dranxo

2

Matrice 400x400! = Grande. De plus, que signifie le plus grand si «il y aura probablement plusieurs valeurs propres de même ampleur»? En terre engourdie: linalg.eig (random.normal (taille = (400 400))) prend environ une demi-seconde. Est-ce trop lent?

meawoppl

@meawoppl oui une demi-seconde est trop lente. C'est parce qu'il fait partie d'un autre algo qui exécute ce calcul plusieurs fois.

puissance

1

@power gotcah. Avez-vous une approximation du vecteur propre. c'est-à-dire qu'elle est probablement similaire à la dernière solution, ou pouvez-vous faire une supposition éclairée sur sa structure?

meawoppl

14

Cela dépend beaucoup de la taille de votre matrice, dans le cas à grande échelle également de sa rareté et de la précision que vous souhaitez atteindre.

Si votre matrice est trop grande pour permettre une factorisation unique et que vous avez besoin d'une grande précision, l'algorithme de Lanczsos est probablement le moyen le plus rapide. Dans le cas non symétrique, l'algorithme Arnoldi est nécessaire, ce qui est numériquement instable, donc une implémentation doit y remédier (est quelque peu difficile à guérir).

Si ce n'est pas le cas dans votre problème, donnez des informations plus précises dans votre question. Ajoutez ensuite un commentaire à cette réponse et je la mettrai à jour.

Edit: [C'était pour l'ancienne version de la question, demandant la plus grande valeur propre.] Comme votre matrice est petite et apparemment dense, je ferais une itération Arnoldi sur B = (IA) ^ {- 1}, en utilisant une initiale la factorisation triangulaire permutée de IA pour avoir une multiplication bon marché par B. (Ou calculer un inverse explicite, mais cela coûte 3 fois plus que la factorisation.) Vous voulez tester si B a une valeur propre négative. En travaillant avec B à la place de A, les valeurs propres négatives sont bien mieux séparées, donc s'il y en a une, vous devriez converger rapidement.

Mais je suis curieux de savoir d'où vient votre problème. Les matrices non symétriques ont généralement des valeurs propres complexes, donc `` plus grandes '' ne sont même pas bien définies. Ainsi, vous devez en savoir plus sur votre problème, ce qui pourrait aider à suggérer comment le résoudre encore plus rapidement et / ou de manière plus fiable.

Edit2: Il est difficile d'obtenir avec Arnoldi un sous-ensemble particulier d'intérêt. Pour obtenir de manière fiable les valeurs propres absolument les plus grandes, vous devez effectuer une itération de sous-espace en utilisant la matrice d'origine, avec une taille de sous-espace correspondant ou dépassant le nombre de valeurs propres censées être proches de 1 ou plus. Sur les petites matrices, cela sera plus lent que l'algorithme QR mais sur les grandes matrices, ce sera beaucoup plus rapide.

Arnold Neumaier
la source

Je dois tester si la plus grande valeur propre est supérieure à 1. La précision ne doit être que de 1e-11. Une matrice "typique" a jusqu'à présent mesuré 386 x 386. J'ai pu y faire une factorisation QR. Cependant, il est également possible d'en avoir beaucoup plus gros. Je commence actuellement à lire sur l'algorithme Arnoldi. Je comprends que cela est lié à Lanczsos.

puissance

Cette information appartient à votre question - veuillez donc la modifier, et ajoutez également plus d'informations (pourquoi les valeurs propres sont-elles réelles? Ou qu'est-ce que la plus grande signifie?) - voir la modification de ma réponse.

Arnold Neumaier

désolé de ne pas m'être expliqué clairement. Je n'ai pas non plus expliqué clairement que les valeurs propres sont complexes. Je teste si des valeurs propres ont une magnitude supérieure ou égale à un.

puissance

1

(I - A)^{- 1}

$(I-A)^{-1}$

1

voir la modification 2 dans ma réponse

Arnold Neumaier

7

$\left|\lambda_{n-1}/\lambda_{n}\right|$

$\lambda_{n-1}$ $\lambda_n$

Pedro
la source

1

et si | λ (n − 1) | = | λ (n) | ?

puissance

@power, l'itération de puissance régulière ne convergera pas. Je ne sais pas dans quelle mesure les méthodes d'extrapolation feront la distinction entre les différentes valeurs propres, vous devrez lire l'article pour cela.

Pedro

2

| λ_{n - 1} | = | λ_{n} |

$|\lambda_{n-1}|=|\lambda_n|$

λ_{n}

$\lambda_n$

λ_{n - 1}

$\lambda_{n-1}$

avez-vous une référence à un article ou un livre académique qui soutient cela? Et si \ lambda_ {n} est complexe?

puissance

5

S'il existe plusieurs valeurs propres différentes de module maximal, l'itération de puissance ne converge que dans des circonstances exceptionnelles. Il oscille généralement d'une manière quelque peu imprévisible.

Arnold Neumaier

5

Il y a eu récemment de bonnes recherches à ce sujet. Les nouvelles approches utilisent des «algorithmes randomisés» qui ne nécessitent que quelques lectures de votre matrice pour obtenir une bonne précision sur les plus grandes valeurs propres. Cela contraste avec les itérations de puissance qui nécessitent plusieurs multiplications matrice-vecteur pour atteindre une grande précision.

Vous pouvez en savoir plus sur la nouvelle recherche ici:

http://math.berkeley.edu/~strain/273.F10/martinsson.tygert.rokhlin.randomized.decomposition.pdf

http://arxiv.org/abs/0909.4061

Ce code le fera pour vous:

http://cims.nyu.edu/~tygert/software.html

https://bitbucket.org/rcompton/pca_hgdp/raw/be45a1d9a7077b60219f7017af0130c7f43d7b52/pca.m

http://code.google.com/p/redsvd/

https://cwiki.apache.org/MAHOUT/stochastic-singular-value-decomposition.html

Si la langue de votre choix n'est pas là, vous pouvez rouler votre propre SVD aléatoire assez facilement; il ne nécessite qu'une multiplication vectorielle matricielle suivie d'un appel à un SVD standard.

dranxo
la source

4

Ici , vous trouverez une introduction algorithmique à l'algorithme Jacobi-Davidson, qui calcule la valeur propre maximale.

Dans cet article, les aspects mathématiques sont explorés. JD permet des matrices générales (réelles ou complexes) et peut être utilisé pour calculer des plages de valeurs propres.

Ici vous pouvez trouver diverses implémentations de bibliothèque JDQR et JDQZ (y compris une interface C, à laquelle vous devriez être capable de lier à partir de R).

Dernier soupir
la source

Je n'ai pu trouver aucune littérature qui déclare explicitement que la méthode Jacobi-Davidson fonctionne pour une matrice réelle et générale.

puissance

À moins que chaque article n'énonce explicitement une restriction et que l'argument de convergence repose sur la restriction qui n'a pas d'importance.

Deathbreath

Voici une autre explication de JD. Les matrices considérées sont complètement générales. Aucune structure spéciale n'est exploitée et les résultats spécifiques aux matrices hermitiennes sont comparés et contrastés, par exemple, la convergence pour les matrices générales est quadratique, mais cubique pour les matrices hermitiennes.

Deathbreath

Merci pour cela. Je ne trouve aucun code C pour une matrice générale, je vais donc devoir écrire le mien. Les liens vers les algorithmes semblent être uniquement pour les matrices hermétiennes.

puissance

1

@power, vous ne trouverez pas non plus dans la littérature un résultat indiquant que les implémentations QR standard convergent pour une matrice réelle et générale - c'est un problème ouvert, et en effet il n'y a pas longtemps, un contre-exemple a été trouvé pour le code QR dans LAPACK.

Federico Poloni

2

Dans votre message d'origine, vous dites:

"J'ai également essayé d'utiliser l'interface ARPACK dans le package igraph R. Cependant, cela a donné une erreur pour l'une de mes matrices."

Je serais intéressé d'en savoir plus sur l'erreur. Si vous pouvez rendre cette matrice accessible au public quelque part, je serais intéressé à essayer ARPACK dessus.

Sur la base de ce que j'ai lu ci-dessus, je m'attendrais à ce qu'ARPACK fasse un très bon travail d'extraction des valeurs propres les plus grandes (ou quelques-unes des plus grandes) d'une matrice clairsemée. Pour être plus précis, je m'attendrais à ce que les méthodes Arnoldi fonctionnent bien dans ce cas et c'est, bien sûr, sur quoi ARPACK est basé.

La convergence lente de la méthode de la puissance lorsqu'il existe des valeurs propres étroitement espacées dans la région d'intérêt a été mentionnée ci-dessus. Arnoldi améliore cela en itérant avec plusieurs vecteurs au lieu de la méthode de l'un en puissance.

Bill Greene
la source

Je vais voir si je peux retrouver mon travail à l'époque. J'ai travaillé là-dessus il y a un an.

puissance

0

Ce n'est pas le moyen le plus rapide , mais un moyen assez rapide consiste à frapper plusieurs fois un vecteur (initialement aléatoire) avec la matrice, puis à normaliser toutes les quelques étapes. Finalement, il convergera vers le plus grand vecteur propre, et le gain de norme pour une seule étape est la valeur propre associée.

Cela fonctionne mieux lorsque la plus grande valeur propre est sensiblement plus grande que toute autre valeur propre. Si une autre valeur propre est d'une amplitude proche de la plus grande, la convergence prendra un certain temps et il peut être difficile de déterminer si elle a convergé.

Dan
la source

1

Merci Dan, cependant: dans mes matrices, certaines des autres valeurs propres auront une magnitude similaire (sinon la même) que la plus grande. Votre méthode est-elle similaire à l'itération de puissance et l'itération de quotient de Rayleigh? Batterson et Smillie (1990) écrivent que pour certaines matrices non symétriques, l'itération de quotient de Rayleigh ne convergera pas. Batterson, S., Smillie, J (1990) "Rayleigh Quotient Iteration for Nonsymmetric Matrices", Mathematics of Computation, vol 55, num 191, P 169 - 178

power

Si d'autres valeurs propres ont la même ampleur que la plus grande ... alors ces valeurs ne sont-elles pas aussi "les plus grandes"?

le

@EMS: Il s'agirait toujours des "valeurs propres les plus importantes", mais la présence de plusieurs valeurs plus importantes tuerait toujours la convergence.

Dan

Je me demande simplement à quelle valeur propre vous voulez qu'il converge. Des choses comme le quotient de Rayleigh / la méthode de puissance sont signifiées quand il y a une plus grande valeur propre distincte. Votre question demande de trouver la plus grande valeur propre, mais il semble que ce ne soit pas vraiment bien défini pour votre problème. Je suis juste induit en erreur par le titre du message.

le

-1

Le package R rARPACK fonctionne pour moi. Et il semble être très rapide car ce n'est qu'une interface pour ARPACK, le package standard pour l'algèbre linéaire clairsemée (ce qui signifie calculer quelques valeurs propres et vecteurs propres).

HoangDT
la source

Bienvenue chez SciComp! Comme l'indique la question, ARPACK ne fonctionne pas pour l'OP, donc cette réponse n'est pas vraiment utile.

Christian Clason

@HoangDT Cette question est antérieure à rARPACK

puissance

Quelle est la façon la plus rapide de calculer la plus grande valeur propre d'une matrice générale?

Réponses: