2017-09-21 22 views
2

Ich habe eine NxN symmetrische und tridiagonale Matrix, die mit einem Python-Code berechnet wurde und ich möchte sie diagonalisieren.Diagonalisierung einer tridiagonalen, symmetrischen spärlichen Matrix mit Python

Im konkreten Fall habe ich es mit N = 6000 zu tun, aber die Matrix kann größer werden. Da es spärlich ist, nahm ich an, dass der beste Weg zur Diagonalisierung darin bestand, den Algorithmus scipy.sparse.linalg.eigsh() zu verwenden, der mit anderen spärlichen und symmetrischen Matrizen (nicht tridiagonal), mit denen ich arbeitete, sehr gut funktionierte. Insbesondere, da ich nur den tief liegenden Teil des Spektrums benötige, gebe ich in der Funktion k=2 und which='SM' an.

jedoch in diesem Fall dieser Algorithmus scheint da nach etwa 20 Minuten nach der Berechnung ich folgende Fehlermeldung nicht zu arbeiten, erhalten:

ArpackNoConvergence: ARPACK error -1: No convergence (60001 iterations, 0/2 eigenvectors converged)

Warum ist das passiert? Ist es ein Problem in Bezug auf einige Eigenschaften von Tridiagonalmatrizen? Welche Python-Routine (und bitte, nur Python!) Kann ich benutzen, um meine Matrix effizient zu diagonalisieren?

Hier ist der angeforderte minimale Code meine Fehler zu reproduzieren:

import scipy.sparse.linalg as sl 
import numpy as np 

dim = 6000 
a = np.empty(dim - 1) 
a.fill(1.) 
diag_up = np.diag(a, 1) 
diag_bot = np.diag(a, -1) 

b = np.empty(dim) 
b.fill(1.) 

mat = np.diag(b) + diag_up + diag_bot 
v, w = sl.eigsh(mat, 2, which = 'SM') 

auf meinem PC die Konstruktion der Matrix 364ms dauert, während der Diagonalisierung der gemeldeten Fehler gibt.

+1

Könnten Sie einen minimalen Arbeits Beispiel nennen? Könnte es ein Problem mit dem Casting geben? 'scipy.sparse' Funktionen arbeiten mit Sparse Arrays, die eine andere Darstellung im Speicher haben und vielleicht hängt das, was Sie beobachten, damit zusammen? Haben Sie versucht, Ihren Code zu profilieren? – norok2

+0

Sie finden [diese] (http://www.sciencedirect.com/science/article/pii/S1063520312001042) sehr nützlich. * O (nlogn) * rockt – Marco

+0

@ norok2 fertig, danke für die Antwort. –

Antwort

2

ARPACK ist gut darin, die großen Eigenwerte zu finden, kann aber Mühe haben, die kleinen zu finden. Glücklicherweise können Sie dies ganz einfach umgehen, indem Sie die in eigsh integrierten Shift-Invert-Optionen verwenden. Siehe beispielsweise here.

import scipy.sparse.linalg as sl 
import scipy.sparse as spr 
import numpy as np 

dim = 6000 
diag = np.empty(dim) 
diag.fill(1.) 

# construct the matrix in sparse format and cast to CSC which is preferred by the shift-invert algorithm 
M = spr.dia_matrix((np.array([diag, diag, diag]), [0,-1, 1]), shape=(dim,dim)).tocsc() 

# Set sigma=0 to find eigenvalues closest to zero, i.e. those with smallest magnitude. 
# Note: under shift-invert the small magnitued eigenvalues in the original problem become the large magnitue eigenvalue 
# so 'which' parameter needs to be 'LM' 
v, w = sl.eigsh(M, 2, sigma=0, which='LM') 
print(v) 

Für dieses spezielle Beispiel Problem, können Sie überprüfen, ob die oben die korrekten Eigenwerte findet, da die Eigenwerte ein explicit formula haben passieren:

from math import sqrt, cos, pi 
eigs = [abs(1-2*cos(i*pi/(1+dim))) for i in range(1, dim+1)] 
print(sorted(eigs)[0:2]) 
Verwandte Themen