Diagonalisation d'une matrice tridimensionnelle symétrique tridimensionnelle avec Python

Question

J'ai une matrice NxN symétrique et tridiagonale calculée par un code Python et je souhaite la diagonaliser.

Dans le cas spécifique, j'ai affaire à N = 6000, mais la matrice peut devenir plus grande. Comme il est clairsemé, j'ai supposé que la meilleure façon de diagonaliser était d'utiliser l'algorithme scipy.sparse.linalg.eigsh(), qui a fonctionné extrêmement bien avec d'autres matrices creuses et symétriques (pas tridiagonales, cependant) avec lesquelles j'ai travaillé. En particulier, puisque je n'ai besoin que de la partie basse du spectre, je spécifie k=2 et which='SM' dans la fonction.

Cependant, dans ce cas, cet algorithme ne semble pas fonctionner, car au bout d'environ 20 minutes de calcul, je reçois l'erreur suivante:

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

Pourquoi est-ce qui se passe? Est-ce un problème lié à certaines propriétés des matrices tridiagonales? Quelle routine Python (et s'il vous plaît, seulement Python!) Est-ce que je peux utiliser pour diagonaliser efficacement ma matrice?

Voici le code minimal demandé de reproduire mon erreur:

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

Sur mon pc la construction de la matrice prend 364ms, alors que la diagonalisation donne l'erreur signalée.

Answer 1

ARPACK est capable de trouver les valeurs propres de grande amplitude mais peut avoir du mal à trouver les plus petites. Heureusement, vous pouvez contourner cela assez facilement en utilisant les options shift-invert intégrées au eigsh. Voir, par exemple, 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) 

Pour ce problème d'exemple particulier, vous pouvez vérifier que ce qui précède est de trouver les valeurs propres correctes puisque les valeurs propres arrive d'avoir un explicit formula:

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