Python scipy eigh (Arpack) gibt falsche Eigenwerte für das verallgemeinerte Eigenwertproblem an

Question

Ich versuche, ein verallgemeinertes Eigenwertproblem mit Arpack zu lösen. Im Moment verwendet der Code LAPACK, aber das ist zu langsam. Wir brauchen nur ein paar Eigenwerte und die Matrizen sind spärlich, daher sollte die Verwendung von Arpack der richtige Weg sein.

Bevor ich mit dem Originalcode arbeite, habe ich beschlossen, einen einfachen Fall mit Scipy Wrapper für Arpack (Eigs) zu testen, aber die Ergebnisse, die ich erhalte, sind falsch und ändern sich jedes Mal, wenn der Code ausgeführt wird.

Mindestarbeitsbeispiel:

    import numpy as np
    from scipy.linalg import eig
    from scipy.sparse.linalg import eigs
    n = 8
    A = np.diag(np.arange(1,n+1,1.0))
    B = np.eye(n) # We want symmetric but a non-diagonal B. eigs gives correct answer for B=np.eye(n)
    B[0][n-1] = 2
    B[n-1][0] = 2


    evals,_ = eigs(A,k=3,M=B,which='LM') 
    print("The eigenvalues obtained by eigs (uses Arpack)")
    print(evals)

    print("Correct eigenvalues using eig (uses Lapack):")
    evals_l,_ = eig(A,b=B)
    print(evals_l)
```

Answer 1

Die Matrix B (M in der Dokumentation) muss gemäß der Dokumentation positiv definieren : "Wenn Sigma None ist, ist M positiv definit" , dies ist zusätzlich zu der ersten Anforderung "M muss eine reelle, symmetrische Matrix darstellen, wenn A ist real " dem dein B folgt. Die Eigenwerte Ihrer aktuellen Matrix B sind -1, 1 und 6. Die Matrix B ist also nicht eindeutig positiv, so dass es naheliegend wäre, dass Eigs nicht funktionieren würden.

Im Vergleich dazu enthält die Dokumentation für eig keine solche Einschränkung für die Matrix B.

Answer 2

Keine wirkliche Antwort auf Ihre Frage, aber die LAPACK-Implementierung von QZ (die das GEP löst) ist bekanntermaßen ziemlich langsam. Mein PR (https://github.com/Reference-LAPACK/lapack/pull/421) behebt das. Sie könnten meine Gabel aus der Quelle bauen und diese verwenden.

Achten Sie auch bei der Verwendung von MKL darauf, dass sie aus irgendeinem Grund den Blockierungsfaktor 1 für die HT-Reduzierung verwenden, was die Leistung erheblich beeinträchtigt. Abhängig von den Einträgen erreicht meine PR eine Beschleunigung von bis zu 250, daher ist es möglicherweise einen Versuch wert.